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I,  INTRODUCTION 

Aerodynamic  design  requirements  demand  the  ability  to  calculate  tables  of 

thermodynamic  properties  and  equilibrium  compositions  of  air  for  thermodynamic 

states  far  from  normal  laboratory  conditions.     Earlier  design  requirements  led 

to  the  production  of  tables  at  very  high  temperatures  ranging  in  density  from 

quite  low  densities  to  intermediate  densities.     The  NBS  tables-  (which  were 

published  under  AEDC  support)  are  examples  of  a  response  to  this  earlier  need. 

These  tables  covered  the  temperature  range  1500  K  <  T  <  15,000  K  and  the 

density  range  10  0  <  c/c     <  10^,  where  p     is  the  density  at  the  standard  condi- 
-        o  o 

tions  of  T  =  273.15  K  and  P  =  1  bar.    More  recently,  there  have  been  indica- 
tions^ of  needs  for  tables  covering  considerably  higher  densities  at  these 
elevated  temperatures . 

In  this  project,  we  were  requested  to  attempt  to  develop  methods  for  the 

calculation  of  the  properties  of  air  for  densities  up  to  1000  times  normal  sea 

level  density.     If  successful,  these  tables,  taken  together  with  the  NBS  tables 

already  published,  would  result  in  tables  being  available  for  the  enormous 
—6  3 

density  range  10      <  c/cq  <  10     (i.e.  nine  orders  of  magnitude)  at  temperatures 
well  outside  those  associated  with  laboratory  experiments.     Such  tables  would 
then  represent  an  extrapolation  from  ordinary  conditions  of  over  two  orders  of 
magnitude  in  temperature  and  at  the  same  time,  an  extrapolation  of  at  least  one 
order  of  magnitude  in  density  from  the  earlier  tables . 

Combined  extrapolations  of  such  magnitudes  in  temperature  and  density  pose 

very  difficult  problems.     They  must  of  necessity  be  based  on  fundamental 

properties  of  the  systems  under  study,  properties  which  might  not  vary  over  the 

range  of  the  extrapolation.     These  properties  must  also  be  used  in  a  framework 

of  fundamental  theoretical  methods  valid  over  the  entire  range  of  temperatures 

and  densities  of  interest.     In  addition,   specific  numerical  methods  must  be 

developed  and  computer  programs  produced  for  transforming  these  numerical  methods 

into  the  actual  calculation  of  tables.     These  numerical  problems  include  the 

often  difficult  task  of  producing  methods  for  the  solution  of  sets  of  non-linear 

3 

algebraic  equations,     with  such  methods  being  required  to  produce  solutions  for 
a  wide  variety  of  values  for  the  unknown  parameters. 
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It  should  be  obvious  that,  because  of  the  extent  of  the  extrapolations 
required,     empirical  correlation  methods  are  not  appropriate  here. 
Such  engineering  methods  involve  the  least  squares  fitting  of  measured  data  to 
equations  containing  arbitrary  parameters.     Such  parameters  generally  do  not 
have  physical  significance  and  can  therefore  only  be  used  with  extreme  caution 
if  at  all  in  extrapolations  beyond  the  range  of  experimental  data  on  which  their 
values  are  based.     In  fact,  using  arbitrary  parameters  obtained  by  least  squares 
fitting,  in  such  long  extrapolations  beyond  the  range  of  their  data  base 
generally  produces  highly  erratic  behavior.     Furthermore,  considerable  addition- 
al difficulty  could  be  anticipated  in  applying  such  methods  to  mixtures.  Mixing 
rules  for  non-physical  parameters  in  the  context  of  empirical  correlations  are 
generally  arbitrary  and  often  not  useful  even  within  the  range  of  the  data  base 
and  hence  their  behavior  in  extrapolations  is  unpredictable. 

The  original  plan  on  which  this  research  was  based  envisioned  the  use  of 

the  then  developing  and  fundamentally  based  integral  equation  theory  for  the 

4 

equation  of  state  of  fluids.      At  the  time  this  work  began,  that  theory,  though 
very  promising  for  one  component  systems,  had  not  been  applied  to  mixtures  nor 
had  any  attempt  been  made  to  integrate  it  into  the  context  of  chemical  or  phase 
equilibrium.     It  was  not  clear  how  one  might  calculate  the  free  energy  and 
chemical  potentials  of  a  multicomponent  mixture  within  this  approach.     The  study 
of  the  possibility  of  using  the  integral  equation  approach  was,  in  fact,  to 
constitute  the  major  part  of  the  research  program  supported  by  this  contract 
with  the  probability  of  success  not  entirely  clear.     A  rather  complex  and 
purely  numerical  approach  was  envisioned  which,  although  it  might  of  necessity 
be  very  complicated  in  a  numerical  sense,  might  nevertheless  be  expected 
ultimately  to  be  made  to  work.     Because  of  its  expected  complexity,  the  basic 
numerical  parts  of  the  problem  were  postponed  in  favor  of  a  close  study  of  the 
details  of  the  integral  equation  approach  itself.     Initially,   this  involved 
examining  integral  equation  methods  as  a  means  for  representing  the  equation  of 
state  of  pure  substances  before  looking  at  them  as  possible  methods  for  mixtures. 

Any  statistical  mechanical  method  which  could  be  used  in  this  work 
(particularly  one  consistent  with  the  earlier  tables)  including  the  integral 
equation  approach,  requires  the  use  of  intermolecular  potential  functions.  For 
the  temperatures  of  interest  in  aerodynamic  problems,  such  intermolecular 
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potential  functions  would  certainly  be  used  beyond  the  range  in  which  data  fits 
were  carried  out  to  obtain  their  paraneters.     This  means  that  methods  used  to 
obtain  these  parameters  would  need  to  be  understood  very  thoroughly.  Further- 
more, because  of  the  implicit  way  in  which  the  potential  functions  appear  in 
statistical  mechanical  theories,  the  effect  of  uncertainties  in  their  parameters 
on  the  accuracy  of  calculations  in  which  they  are  used  is  not  always  clear.  For 
these  reasons,  much  of  the  time  in  this  research  program  had  to  be  spent  in  the 
study  of  the  means  by  which  intermolecular  potential  functions  are  obtained 
from  experimental  data  of  various  kinds  and  in  the  study  of  the  effect  on  cal- 
culated tables  of  extrapolating  incorrect  intermolecular  functions  .  Although 
this  extrapolation  of  the  potential  parameters  beyond  the  range  of  the  data 
could  be  expected  to  introduce  errors,  these  errors  could  be  expected  to  be  far 
less  than  those  obtained  from  extrapolation  of  least  squared  fits  to  data  such 
as  are  common  in  the  usual  empirical  approaches.        Furthermore,  one  might  expect 
the  general  behavior  of  the  functions  calculated  to  be  reasonable.     This  follows 
from  the  fact  that  the  essential  behavior  of  the  potential  function  would  be 
correct  and  from  the  fact  that  the  "actual"  potential  functions  would  not  vary 
with  temperature  so  that  the  relationship  between  the  potentials  chosen  and  the 
"real"  ones  would  not  contain  any  hidden  surprises  on  extrapolation  to  higher 
temperatures . 

Although  the  present  work  emphasizes  the  high  density  region,  it  had  to  be 
designed  to  retain  the  earlier  tables.     Since  these  earlier  calculations  already 
covered  densities  up  to  100  times  normal  sea  level  density,  the  extension  from 
100  to  1000  times  normal  density  needed  to  be  done  using  methods  which  included 
the  earlier  approach.     At  the  highest  densities,  the  earlier  tables  included  a 
correction  for  the  second  virial  coefficient  based  on  the  intermolecular  forces 
between  the  molecules.     Hence,  almost  any  approach  is  consistent  with  the 
earlier  tables  if  it  makes  use  of  a  statistical  mechanical  theory  based  on  these 
same  intermolecular  forces  provided  that  this  theory  has  the  same  (and  correct) 
linear  term  in  its  low  density  expansion. 

As  the  work  progressed,  alternative  approaches  were  examined  in  an  attempt 
to  avoid  the  numerical  complexity  expected  in  the  integral  equation  approach. 
These  other  methods  were  rejected  either  because  their  complexity  did  not  offer 
any  advantage  over  the  integral  equation  approach  or,  more  usually,  because  they 
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required  the  extrapolation  of  parameters  which  had  no  physical  basis.  One 
alternative  which  appeared  to  hold  promise,  resulted  from  research  by  one  of  us 
(Lester  Haar) .     An  equation  of  state  for  single  component  fluids  was  developed 
which  seemed  to  us  capable  of  being  used  in  this  work  in  a  computational 
procedure  far  simpler  than  that  projected  for  the  integral  equation  method. 
Also,   it  was  based  on  fundamental  ideas  and  was  thought  to  provide  a  natural 
extension  of  the  previous  NBS  tables.     The  approach  is  based  on  an  equation  of 
state  due  to  Kaar  and  Shenker  ,    (HS) ,  which  was  developed  for  the  extrapolation 
of  low  density  experimental  data  of  one  component  system  to  high  densities  along 
isotherms  at  ordinary  temperatures.     According  to  the  assumptions  on  which  it  was 
based,   the  HS  equation  was  expected  to  work  well  at  temperatures  above  the 
critical  temperature  and  to  improve  with  increasing  temperature  and  this  was 
indeed  found  to  be  the  case  by  them.     They  also  found  indications  that  the  method 
could  be  used  to  produce  engineering  calculations  of  reasonable  accuracy  even 
below  the  critical  temperature  into  the  liquid  range. 

The  HS  equation  uses  second  virial  coefficient  data  at  each  temperature  as 
a  basis  for  extrapolating  to  high  densities.     In  our  studies  of  the  relationship 
between  second  virial  coefficients  and  intermolecular  potentials,  we  saw  that 
the  determination  of  intermolecular  forces  from  second  virial  coefficients,  if 
done  properly,  could  be  used  to  produce  a  temperature  extrapolation  of  the 
second  virial  coeffcient  to  higher  temperatures  even  as  high  as  required  in 
aerodynamic  calculations.     Thus,  given  an  intermolecular  potential  function  for 
a  particular  substance,  a  complete  high  temperature  PVT  surface  could,  in 
principle,  be  produced  for  that  substance,  with  the  potential  function  being 
used  to  calculate  second  virial  coefficients  at  all  temperatures  of  interest 
and  these  virial  coefficients,  in  turn,  serving  as  the  basis  for  the  HS  equation 
at  all  densities. 

Before  it  could  be  used  in  this  work,  the  HS  equation  (which  was  developed 
for  a  pure  substance)  had  to  be  extended  to  mixtures  in  the  context  of  a 
reacting  gas.     This  we  have  done  and  the  details  of  the  formalism  will  be 
presented  in  this  report.     The  resulting  formalism  was  used  by  us  to  produce 
the  attached  set  of  tables.     A  difficulty  is  expected  to  arise  in  this  formalism 
whenever  there  is  a  need  for  intermolecular  potential  functions  between 
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fragments  which  are  present  in  the  mixture  at  high  temperatures  but  which  do  not 
exist  at  ordinary  temperatures.     Potential  functions  between  such  molecular 
fragments  can,  however,  often  be  characterized  approximately  by  means  of  the 
potential  function  which  describes  the  vibrational  spectra  of  the  combined 
fragments  or,  failing  such  spectra,      from  estimates  of  molecular  sizes.  Since 
estimates  of  the  intermolecular  forces  between  the  various  species  had  to  be 
made  for  the  earlier  tables   (including  those  species  not  present  at  ordinary 
temperatures),  the  basic  data  needed  for  this  approach  were  actually  already 
available  at  the  start  of  this  work.     In  order  to  ensure  consistency  in  the 
initial  calculations,  we  decided  to  use  these  same  estimates  of  these  forces  in 
the  calculations  on  which  we  report  here. 

Although  the  formal  framework  of  the  approach  used  by  us  is  a  valid  one,  the 
values  contained  in  the  attached  tables  should,  for  a  number  of  reasons,  only  be 
considered  an  interim  set  of  values :- 

1.  The  estimates  of  the  intermolecular  potentials  used  are  based  on  an 
old  analysis  by  Wbolley    in  terms  of  the  (12,6)  potential  function.     As  a 
result  of  our  earlier  work  under  AEDC  support,  better  potentials  are  now  known 
and  should  be  used  for  the  major  species  and  these  same  potentials  should  also 
be  used  in  place  of  the  (12,6)   for  the  new  estimates  of  those  associated  with 
the  minor  species.     As  already  mentioned,  we  used  these  older  estimates  so  as 
to  ensure  consistency  between  these  tables  and  the  earlier  tables. 

2.  In  several  cases,  certain  unknown  intermolecular  potential  functions 
were  arbitrarily  taken  as  equivalent  to  others  that  were  known.  The  effect  of 
this  could  be  important  at  the  very  highest  densities. 

3.  As  we  shall  describe  below,  we  have  used  an  ad  hoc  approach  in  the 
determination  of  hard  sphere  diameters  at  those  temperatures  for  which  the 
slope  of  the  second  virial  coefficient  is  negative. 

4.  In  this  work  only  the  compressibility  factor  and  equilibrium  composi- 
tions have  been  calculated. 

Each  of  these  compromises  has  been  made  in  order  to  produce  results  to  test  the 
approach,  leaving  the  production  of  more  extensive  tables  for  possible  later  work. 
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Grabau  and  Brahinsky  ,   (GB) ,  attempted  an  extrapolation  of  the  earlier  NBS 

3 

tables  to  10     times  normal  density  for  air  and  nitrogen.     Their  method  was 
essentially  semi-empirical.     Their  method  was  based  on  an  equation  having,  in 
part,  a  fundamental  basis  and  which  contained  a  dependence  on  ideas  partially 
related  to  the  use  of  intermolecular  forces.     Their  method  could  not  be  applied 
to  a  gas  of  varying  composition  (i.e.  a  reacting  gas),  for  which  reason  they 
restricted  their  air  calculations  to  temperatues  below  6000K.     Their  method 
would  also  be  expected  to  lose  accuracy  very  rapidly  at  any  temperature  as  the 
density  increased  due  to  a  need  for  successive  subtractions  and  because  of  their 
neglect  of  higher  order  terms.     Since  their  method  was  based  mainly  on  a  graph- 
ical extrapolation,  even  qualitative  estimates  of  its  adequacy  could  not  be  made 
by  them. 

In  this  report,  we  compare  our  results  with  the  extrapolation  of  Grabau 
and  Brahinsky.     These  comparisons  are  interesting  and  will  be  discussed  below. 
The  agreement  between  the  two  extrapolations  at  intermediate  density  at  the 
temperatures  which  they  have  in  common  is  expected  since  the  GB  extrapolations 
were  based  on  the  earlier  NBS  tables  with  which  our  results  are  entirely 
consistent.     The  two  results  show  considerable  disagreement  at  the  highest 
densities  for  reasons  discussed  below. 

Our  report  contains  plots  at  two  representative  temperatures  of  species 
concentrations  for  some  important  species.     These  plots  show  some  interesting 
and  possibly  unexpected  density  behavior  as  discussed  below.     Information  on 
composition  behavior  is  not  available  in  the  approach  of  Grabau  and  Brahinsky 
since  their  method  is  applied  to  the  thermodynamic  properties  only.     Of  partic- 
ular interest  in  our  plots  are  the  dependences  of  concentrations  on  density  as 
produced  in  three  widely  used  approximations  -  the  ideal  gas,  the  second  virial 
coefficient  gas  and  the  HS  gas.     Regardless  of  the  ultimate  accuracy  of  our  own 
model,   the  differences  among  these  three  approximations  can  be  expected  to  be 
indicative  of  the  magnitude  of  the  effect  of  the  analogous  three  approximations 
in  any  other  (and  possibly  more  accurate)  theoretical  model. 

Many  unexpected  technical  numerical  and  computer  programming  problems  were 
encountered  in  adding  on  the  extreme  high  density  end  of  the  calculation. 
These  had  to  do  with  numerical  difficulties  in  obtaining  solutions  and  will  be 
discussed  only  in  passing  in  this  report. 
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II.     THE  EQUATION  OF  STATE  OF  HAAR  AND  SHENKER 

A.  BACKGROUND 

The  first  step  in  the  calculation  of  the  thermodynamic  properties  of  a  re- 
acting mixture  involves  the  computation  of  the  composition  corresponding  to 
reaction  equilibrium.  Within  the  framework  of  the  formalism  previously  used  by 
us,  this  corresponds  to  the  solution  of  the  equations  associated  with  the  law  of 
mass  action,  as  modified  to  take  account  of  any  non-ideal  effects,  subject  to 
the  conservation  laws  for  nuclear  types.     The  mass  action  relations  can  be 
written 

C.     =    K.(p/porWi  n  V.WJl  C,Vik  (1) 

1       1  I  1    k  k 


where  the  C,    are  the  concentrations  of  the  reference  species,  v.,  the 

K-  ~  co^  ik 

stoichiometric  coefficients  for  the  reaction,  K.  =  K.(T/T  )      with  K.  the 

'     1         l         o  1 

equilibrium  constant,  -co^  =  Ev^-l  the  net  production  of  particles  across  the 

reaction,  p     the  density  at  standard  conditions  (for  P  =  one  atmosphere 

°  (O 
pressure  and  for  T  =  T    =  273.15  kelvins)  and  where  y*       is  the  effective 

°  th  th 

activity  coefficient  for  the  I      non-ideal  effect  for  the  i      reaction.  A 

detailed  discussion  of  this  equation  and  its  derivation  are  contained  in 

reference  (3)  and  will  not  be  repeated  here.     The  y^  associated  with  the  Debye- 

Huckel  theory  and  that  associated  with  the  second  virial  coefficient  are 

contained  in  Appendix  B  of  reference  (3) .     The  derivation  of  the  y^  for  the  HS 

equation  constitutes  a  major  part  of  this  report.     For  the  moment,  it  will  be 

enough  to  state  that,  in  principle,  a  y!^  can  be  obtained  for  any  equation  of 

state  but,  in  practice,  the  procedure  is  very  complicated  and  not  always  clear. 

The  understanding  that  formal  expressions  for  the  y!  for  the  HS  equation  can  be 

obtained  constitutes  the  motivation  for  the  following  extended  discussion  of 

that  equation  of  state,  of  its  development  and  of  tests  of  its  validity. 
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Haar  and  Shenker  (HS) ,  developed  an  extremely  simple  equation  of  state 
based  on  the  virial  expansion  and  on  the  behavior  of  the  virial  coefficients  at 
high  temperatures.     The  HS  equation  requires  only  a  knowledge  of  the  second 
virial  coefficient  and  its  first  derivative  at  each  temperature  and  from  this 
the  behavior  at  all  densities  is  obtained.     As  already  mentioned,  the  HS 
equation  can  be  used,  in  a  very  simple  and  straightforward  manner,  to  extra- 
polate PVT  data  in  both  the  temperature  and  density  directions.     It  is  also 
expected  to  improve  with  increasing  temperatures  since  the  validity  of  the 
assumptions  on  which  it  is  based  improves  with  increasing  temperature.     It  is 
therefore  particularly  well  suited  for  these  calculations. 

The  development  of  the  equation  of  state  was  motivated  by  earlier  work  by 
9 

Haar  and  Levelt  Sengers    who  showed  that  only  two  parameters  are  required  to 
correlate  thermodynamic  properties  for  a  number  of  simple  non-polar  gases  along 
individual  isotherms.     The  two  parameters  of  Haar  and  Levelt  Sengers  had  to  be 
different  for  each  temperature  and  were  obtained  by  fits  to  experimental  data  on 
each  isotherm.     These  two  parameters  could,  in  principle,  be  determined  from 
experimental  data  at  each  temperature  in  a  number  of  different  ways. 

The  equation  of  state  of  Haar  and  Shenker  also  has  two  parameters  for  each 
isotherm  with  the  values  of  the  parameters  being  determined  from  values  of  the 
experimental  second  virial  coefficient.     Since  there  are  two  parameters  at  each 
temperature,  two  properties  of  the  second  virial  coefficient  are  required  at 
each  temperature  to  determine  their  values.     Haar  and  Shenker  chose  for  these 
the  second  virial  coefficient  and  its  first  derivative.     The  parameters  so 
obtained  can  be  regarded  as  being  an  effective  temperature  dependent  molecular 
size,  which  sets  the  scale  of  density  on  the  isotherm,  and  an  effective  temp- 
erature dependent  molecular  well  depth  which  sets  the  temperature  scale  on  which 
the  isotherms  are  assigned.     Each  of  these  parameters  is  associated  with  the 
experimental  system.     In  the  following  derivation  these  two  quantities  will  be 
considered  to  be  only  slowly  varying  functions  of  the  temperature  which  greatly 
simplifies  their  calculation.     This  assumption  is  certainly  valid  in  the  range 
of  temperatures  of  interest  in  this  work.     In  this  way  Haar  and  Shenker 
developed  an  equation  of  state  for  the  correlation  and  prediction  of  high 
density  data  at  temperatures  above  critical  using  two  temperature  dependent 
parameters.     The  parameters  chosen  have  a  fundamental  basis  and  a  simple  method 
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for  their  determination  was  devised.     The  fact  that  the  parameters  so  chosen 
could  also  be  associated  with  notions  of  corresponding  states  promised  general 
success . 

The  repulsive  energy  between  two  molecules  in  a  gas  can  generally  be 

characterized  by  a  very  steep  function  at  small  intermolecular  separations.  To 

a  first  approximation  this  repulsion  can  be  described  by  the  interaction  between 

hard  spheres.     Since  repulsive  effects  are  known  to  dominate  at  high  temper- 
10 

atures     ,  it  is  reasonable  to  take  the  equation  of  state  of  a  gas  of  hard  spheres 
(for  which  there  are  only  repulsive  effects)  as  a  starting  point  in  the  develop- 
ment  of  an  equation  of  soaoe  for  any  gas  at  hit-,  temperatures.     To  facilitate 
this,  Haar  and  Shenker  express  the  actual  equation  of  state  as  the  sum  of  a 
hard  sphere  contribution  (to  be  calculated  by  a  method  as  yet  unspecified)  plus 
the  difference  between  the  actual  equation  of  state  and  this  hard  sphere 
contribution.     This  involves  no  approximations  since  the  two  parts  add  up 
identically  to  the  actual  equation  of  state  regardless  of  how  the  hard  sphere 
contribution  is  handled.     The  first  approximating  assumption  consists  in  taking, 
for  the  hard  sphere  contribution,  the  result  obtained  from  Percus-Yevick  theory^ 
using  the  compressibility  equation  of  state.     It  has  been  shown  that  this 

representation  differs  only  slightly  from  exact  hard  sphere  theory  up  to 

12 

densities  approaching  2/3  the  close  packing  density     .     A  second  (and  more 
serious)  approximation  involves  the  choice  of  method  for  obtaining  the  hard 
sphere  diameters  needed  in  the  Percus-Yevick  theory. 


The  derivation  of  the  HS  equation  of  state  starts  from  the  assumption  that 

the  N  body  potential  of  the  fluid  can  be  represented  by  a  sum  of  pair-potentials 

13 

and  that  the  Ursell-Mayer  virial  expansion  in  the  density      is  valid  for  all 
potential  functions  of  interest.     The  virial  series  for  any  potential  is  then 
transformed  into  a  rapidly  convergent  expansion  about  the  hard  sphere  series. 
Finally,  the  equation's  parameters  are  fixed  by  imposing  as  boundary  conditions 
the  requirement  that  the  first  correction  to  the  ideal  gas  be  valid.     This  last 
follows  automatically  when  the  second  virial  coefficient  is  used  to  determine 
the  parameters  and  is  the  basis  of  the  consistency  between  our  model  and  the 
earlier  NBS  tables. 
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B.       DERIVATION  FOR  A  PURE  SUBSTANCE 

The  Ursell-Mayer  expansion  for  the  equation  of  state  is  written, 


n=l 

PS 

where  Z  is  the  compressibility  factor,  defined  by  Z  =  —  ;  p  the  reciprocal 
volume,  P  the  pressure,   8  =  1/kT,  and  N  the  number  of  molecules  in  the  system. 
The  density  expansion  (2)   for  the  compressibility  factor  is  certainly  valid  in 
the  gas  phase  and  can  be  considered  to  be  an  exact  representation  of  the  com- 
pressibility factor  Z.     The  B     in  (2)  are  the  virial  coefficients  and  are  well- 

14  n 
known      integrals,  obtained  from  statistical  mechanics,  involving  the  inter- 
molecular  interactions  among  2,3,4,  etc.  molecules  respectively.     We  now 
formally  develop  each  of  these  virials  about  that  for  a  hard  sphere  of  some 
Cas  yet  arbitrary)  diameter  via  the  identity 


B    =  Bh-S'  +  B    -  Bh-S>         ,  (3) 
n        n  n  n 


where  B^*S*   is  the  virial  coefficient  for  the  hard  sphere.     Using  (3),  we 

can  rewrite  eq.  (2) 


i       "  -,        n  n 

n=l  n=l 


where,  consistent  with  the  formulation,  B^"       =  B^  =  1.     It  should  be  noted  that 
(4)  is  an  identity  and  so  does  not  involve  any  new  assumptions.     (4)  is  there- 
fore still  an  exact  representation  for  Z.     The  first  sum  on  the  right-hand  side 
of  Eq.    (4)  is  the  equation  of  state  that  would  be  obtained  for  a  gas  made  up  of 
identical  hard  sphere  molecules.     The  hard-sphere  gas  has  been  studied  exten- 
sively in  computer  "experiments"  via  molecular  dynamics"'""'  and  Monte  Carlo"^  cal- 
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dilations  and  theoretically  via  the  Percus-Yevick  approximation .        The  hard- 
sphere  summation  in  (4)  can  therefore  be  considered  to  be  known  for  any  partic- 
ular sphere  size.     The  task  inherent  in  the  evaluation  of  (4)  is  to  obtain  a 
simple  representation  for  the  perturbation  terms,  i.e.  the  second  sum  on  the 
right  of  Eq.    (4),  and  to  obtain  a  proper  hard  sphere  diameter  to  use  in  the 
first  term. 

We  now  present  an  argument  which  shows  that  the  perturbation  terms,  i.e. 
the  second  sum  on  the  right-hand  side  of  (4) ,  converge  rapidly  above  the  gas- 
liquid  critical  temperature  so  that  at  such  temperatures  only  the  term  linear  in 

density  needs  to  be  considered.     The  prospect  that  this  convergence  might 

17 

persist  to  somewhat  lower  temperatures  is  implied  in  work  by  Woolley 

To  illustrate  our  argument,  in  Fig.  1  we  plot  a  few  of  the  lower  virial 

coefficients  for  the  Lennard- Jones  0-2,6.)  pair-potential.     The  reduced  virials 

B„,  B„,  and  B.  with 
2      3  4 

*  ..2      „  3,n-l 

B      =B  /  [-r  ttNo  J  , 
n         n  3 


are  plotted  against  the  reduced  temperature  T    =  (3e)     ,  where  a  is  the  Lennard- 
Jones  length  parameter  and  e  the  well  depth.     At  low  temperatures  the  contribu- 
tion from  the  attractive  part  of  the  potential  (regions  of  negative  energy)  is 
important  and  the  lower  virials  tend  to  large  negative  values.     At  the  higher 
temperatures,  the  repulsive  part  of  the  potential  (regions  of  positive  energy) 
tends  to  yield  the  dominant  contribution  and  the  virials  become  positive.  Thus, 
in  Fig.   1  we  see  that  B^  monotonically  increases,  with  decreasing  slope,  from 
large  negative  values  to  positive  ones,   finally  passing  through  a  maximum.  It 
is  apparent  that  as  the  temperature  increases,  the  contribution  of  attraction 
decreases  relative  to  that  of  repulsion,  even  though  the  contribution  of 
repulsion  is  itself  slowly  decreasing  with  increasing  temperature.     We  suggest 

that  at  some  temperature,  not  much  beyond  the  temperature  where  B    achieves  its 

*  * 
maximum,  the  contribution  of  attraction  to  B    is  relatively  small.     This  occurs 

*  r 

at  a  reduced  temperature  somewhat  above  T    =30.     Thus  if  the  hard-sphere 
contribution  is  chosen  appropriately,  we  would  expect  that  the  perturbation 
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terms  linear  in  density  would  be  quite  small  at  temperatures  somewhat  above 
T    =30.     Likewise  we  expect  that  for  the  third,  fourth,  etc.  virials  the 
contribution  of  attraction  becomes  relatively  small  at  temperatures  starting 
somewhat  above  the  temperatures  for  which  these  virial  coefficients  achieve 
their  (initial)  maxima.     The  maximum  for  the  third  virial  coefficient  occurs  at 


about  T    =1.25  corresponding  to  a  temperature  near  the  liquid-vapor  critical 

point  for  simple  substances.     The  temperature  at  which  the  fourth  virial 

coefficient  would  achieve  its  maximum  is  also  in  this  neighborhood.     (This  is 

18 

also  true  for  the  fifth  virial  coefficient  not  shown  in  the  figure     ).We  invoke 
a  corresponding  states  argument  and  assume  that  the  disappearance  of  attractive 
contributions  to  the  third  and  higher  virial  coefficients  above  the  critical 
temperature  should  be  a  general  property  of  any  simple  gas.     Thus,  we  suggest 
that  at  temperatures  somewhat  above  the  critical  temperature  the  second  virial 
coefficient  includes  all  of  the  major  effects  due  to  molecular  attractions,  and 
that  the  higher  virial  coefficients  are  primarily  determined  by  the  repulsive 
interaction.     For  each  substance,  at  any  given  temperature,  we  therefore  repre- 
sent the  repulsion  between  the  molecules  by  that  in  a  gas  made  up  of  identical 
hard  spheres  whose  diameter  is  somehow  chosen  so  as  to  be  appropriate  to  that 
substance  at  that  temperature.     We  then  use  this  hard  sphere  gas  to  represent 
the  total  contribution  of  each  of  the  virials  above  the  second  at  that  tempera- 
ture.    This  certainly  should  be  a  good  model  at  the  temperatures  of  interest  in 


Based  on  the  preceding  arguments,  for  temperatures  of  interest  here,  the 
second  sum  in  the  equation  of  state,  Eq.    (4),  can  be  truncated  after  the  term 
linear  in  the  attractive  contribution.     The  hard  sphere  part,  on  the  other  hand, 
contains  contributions  to  all  orders  of  the  density.     To  represent  this  hard 
sphere  part,  we  employ  the  results  of  the  Percus-Yevick  theory"^  using  the 
compressibility  form  for  the  equation  of  state.     This  is  a  good  approximation 
at  low  densities,  and  is  in  error  by,  at  most,  5%  at  densities  approaching  2/3 
close  packing  of  hard  spheres,  a  density  well  beyond  those  of  interest  here. 
The  equation  of  state  is  written,  therefore, 


* 


our  work. 


2 


+  4y   (^  -  1)  , 


(5) 
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where  4y  =  bp,  and  b  is  the  hard  sphere  second  virial  coefficient, 


tt  N  a 


3 


(6) 


a  being  the  temperature  dependent  hard  sphere  diameter.     For  a  given  temperature 
Eq.   (5)  is  a  two  constant  equation  of  state,  these  constants  being  the  hard 
sphere  diameter  a  and  the  well  depth  associated  with  the  representation  of  . 
An  important  feature  to  be  used  below  is  that  the  equation  (5)  is  easily  inte- 
grable  in  closed  form  to  yield  a  free  energy. 

As  already  mentioned,  to  evaluate  the  two  equation  of  state  parameters  for 
a  particular  gas  at  each  temperature  we  shall  employ  the  numerical  values,  for 
that  gas,  of  the  second  virial  coefficient  and  its  first  temperature  derivative. 
To  bring  out  the  connection  between  these  two  parameters  and  the  molecular 
diameter  and  the  intermolecular  well  depth,  we  introduce  an  effective  inter- 
molecular  interaction  that  has  the  general  features  of  a  typical  pair  potential, 
except  that  it  is  specifically  characterized  by  a  hard  sphere  cut-off  at  some 
diameter  a(T)  at  which  point  it  is  joined  to  an  attractive  bowl  of  well  depth 
e(T)  CNote  that  we  have  explicitly  indicated  the  temperature  dependence  of  these 
quantities.)     The  purpose  of  this  effective  function  is  to  provide  a  means  for 
transforming  the  repulsive  and  attractive  parts  of  the  "actual"  potential  of 
the  gas  into  an  explicit  hard  sphere  diameter  and  a  well  depth.  Typical 
functions  of  this  kind  are  shown  in  Fig.  2. 

A  simple  numerical  method  has  been  developed  for  the  determination  of  the 
parameters  a  and  e  at  each  temperature.     The  method  starts  with  tables  of 


the  effective  potential  function  (i.e.  the  function  with  the  hard  sphere  cut- 
off).    At  any  given  temperature,  these  reduced  quantities  are  required  to  yield 
the  values  associated  with  the  experimental  system  being  described.     Thus  the 
conditions  on  these  reduced  quantities  are, 


reduced  second  virial  coefficients  B    and  their  first  derivative  T    — ,  for 


dT 


B 


exp 


(T)  =  b  B  (kT/e) 


(7) 
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and  dBexD^  *  dB*  (kT/e) 

T  —^2-        =  b  T    2£_  (8) 

dT 

where  T  is  the  temperature  of  interest.     On  the  right-hand  side,  we  have 
indicated  the  dependence  on  T    by  kT/e  to  emphasize  the  fact  that  e/k  is  an 
unknown  quantity.     The  right-hand  sides  in  (7)  and  (8)  are  the  reduced  quanti- 
ties as  calculated  for  the  effective  potential  function  (i.e.  the  potential  with 

the  hard  core).     It  should  be  noted  that  the  right  hand  side  of  equation  (8)  contains 

da 

the  implicit  assumptions  that  —  and  de/dT  can  be  neglected  in  the  calculation  of 
a  and  e.     On  dividing  Eq.   (8)  by  Eq.   (7),  there  results 

T         ^exp_         T*      dB*.  _  0* 

B  dT  *  *  (9) 

exp  B  dT 


* 

It  is  a  simple  matter  to  produce  a  table    of  values  for  the  quantity  Q    as  a 

*  * 
function  of  T    for  the  effective  potential  function.     By  way  of  illustration  Q 

values  are  listed  in  Table  1  for  a  particular  effective  potential.     The  same 

quantity  (i.e.  the  left  hand  side  of  (9))  is  then  calculated  from  experimental 

data  as  a  function  of  T.     (9)  is  then  solved  for  e/k  at  a  given  experimental 

* 

temperature  T.     This  is  done  by  starting  with  the  experimental  value  of  Q  at 

* 

that  value  of  T,  and,  by  interpolation,  finding  that  value  of  Q    and  the  value 
of  T    associated  with  it  in  the  table  of  values  calculated  for  the  effective 
potential.     e/k  is  then  calculated  at  that  temperature  from  the  relation 


e/k  =  T/T* 

This  procedure  guarantees  that  Eq.    (9)   is  automatically  satisfied  at  the  chosen 

experimental  value  of  T.     This  e/k  value  is  then  used  in  (7)  to  obtain  b  and 

3 

therefore  a  .     By  carrying  out  this  procedure  at  each  temperature,  a(T)  and  e(T) 
are  determined  for  all  experimental  points.     The  equation  of  state  is  then  cal- 
culated  at  each  temperature  T  using  b(T)  and  B     (T  )  in  equation  (5).     (The  use 


of  B^  (T  )  is  equivalent  to  the  use  of  e(T)) 
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One  of  the  difficulties  associated  with  the  above  procedure  involves  the 

calculation  of  dB/dT  from  experimental  data  which  are  generally  not  smooth  and 

not  presented  at  convenient  temperature  intervals.     A  reasonable  way  to  do  this 

is  first  to  fit  the  experimental  B(T)   data  to  a  realistic  (i.e.  one  which  contains 

a  dependence  on  r  for  small  r)  potential  function  (as  opposed  to  the  effective 

potential  which  has  a  hard  core  cut-off  for  small  r)  and  to  calculate  smooth  tables 
exp 

of  B        and  — — values  using  that  function, 
exp  dT  6 

The  procedure  outlined  for  solving  Eq.    (9)  does  not  work  for  those  temp- 
erature at  which  dB      /dT  <  0.     Because  the  effective  potential  has  a  hard  core 
exp  r 

repulsion,  its  second  virial  coefficient  does  not  have  a  region  of  negative 
slope  and  solution  of   (9)  becomes  impossible  since  the  negative  value  of  Q 
associated  with  the  experimental  system  is  being  sought  in  the  table  for  the 
effective  potential  which  contains  only  positive  values.     Under  such  conditions, 
we  have  proceeded  by  neglecting  the  attractive  contribution  by  setting  B  equal 
to  unity  in  (7).     This  leads  to  a  negligible  discontinuity  in  b(T)  but  in  a  non- 
negligible  one  in    its  temperature  derivative. 

C.       COMPARISONS  WITH  EXPERIMENTAL  DATA  AND  OTHER  THEORIES 


1.       Sensitivity  to  the  selection  of  the  attractive  part  of  the  effective 
potential. 

It  should  be  obvious  that  the  effective  potential  plays  no  fundamental  role 
in  the  HS  theory  but  is  used  only  for  computational  convenience,  being  used  to 
extract  an  effective  hard  sphere  diameter  from  the  second  virial  coefficient. 
It  clearly  should  not  be  allowed  to  introduce  any  of  its  own  character  into  the 
calculation.     For  this  reason,  we  shall  precede  a  detailed  comparison  with 
experiment,  by  an  examination  of  the  sensitivity  of  our  method  to  any  particular 
choice  for  the  shape  of  the  bowl  used  in  the  effective  potential.     In  Fig.   3  a 

■k 

plot  of  the  sphere  diameter  a(T)  vs.  T    is  shown  for  the  two  effective  potential 
functions  corresponding  to  m=9  and  m=10 .     As  expected  from  the  behavior  of  a 
typical  intermolecular  potential  function  at  small  separation,   the  effective 
sphere  diameter  is  a  monotonically  decreasing  function  of  the  temperature.  The 
scale  of  the  abscissa  is  normalized  so  that  a(T)  =  1  for  the  effective  potential 
function  given  by  m=12.     The  parameters  a  and  z  were  calculated  from  equations 
(7)  and  (9)   for  each  of  these  effective  potentials  using  values  for  B„  and  dB/dT 
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calculated  with  a  Lennard-Jones  (12,6)  potential  in  place  of  experimental  data. 
The  choice  of  a  particular  effective  potential  function,  that  is,  the  choice  for 
the  shape  of  the  potential  bowl  in  Figure  2,  is  far  from  unique,  for,  in 
addition  to  the  dependence  on  m,  the  bowl  shape  could  be  affected  by  adjusting 
the  exponent  in  the  attractive  term. 

At  a  reduced  temperature  T    =3.0,  a(T),  and  hence  the  equation  of  state, 
was  found  to  be  relatively  insensitive  to  m,  that  is,  to  the  shape  of  the  bowl 
appended  to  the  hard  sphere  core.     At  higher  temperatures,  the  sphere  diameter 
becomes  completely  independent  of  m.     In  such  a  case,  the  procedure  used  for  the 
evaluation  of  the  sphere  size  can  be  further  simplified,  as  will  be  discussed 
below.     On  the  other  hand,  as  the  temperature  is  reduced,  the  hard-sphere  size 
tends  to  become  increasingly  sensitive  to  m.     This  results  in  a  useful  procedure 
for  determining  an  optimum  value  for  m  for  a  given  substance.     This  consists  of 
comparing  experimental  PVT  data  with  those  predicted  by  the  HS  equation  at  a 
low  temperature  for  effective  potentials  characterized  by  several  values  of  m 
until  a  best  fit  is  obtained.     Because  of  the  insensitivity  to  m  already 
described,     the  value  of  m  chosen  can  obviously  be  used  at  higher  temperatures. 
The  method  has  thus  been  modified  to  produce  a  good  low  temperature  fit.  Haar 
and  Shenker  used  an  isotherm  near  3/2  times  the  critical  temperature  for  the  low 
temperature  fitting. 

We  include,  in  the  next  section,  a  comparison  of  the  predictions  of  the  HS 
equation  with  other  theories  and  with  experimental • data  for  argon  and  nitrogen. 
It  is  easily  determined  that  m=9  and  m=10,  respectively,  are  reasonable  values 
of  m  for  these  fluids.     Since  the  sensitivity  of  the  equation  of  state 
properties  to  m  is  weak  except  at  low  temperatures,  m=9  can  also  be  taken  for 
nitrogen  when  the  temperatures  of  interest  do  not  extend  much  below  twice  the 
critical. 

2 .       Relationship  to  other  fundamental  equations  of  state. 

It  is  useful  to  compare  the  HS  equation  with  equations  of  state  which  have 
the  appearance  of  being  more  fundamentally  based.     The  latter  are  invariably 
much  more  complicated  than  is  the  HS  equation  so  would  have  to  produce  far 
superior  results  to  it  to  warrant  their  choice  over  it.     In  this  comparison, 
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we  shall  include  an  example  of  an  integral  equation  (the  Percus-Yevick)  as  well 
as  the  perturbation  theory  approach  of  Barker  and  Henderson.     Although  the  latter 
is  developed  totally  within  the  language  of  statistical  mechanics,  it  is  not 
unrelated  to  the  Haar-Shenker  approach  which  has  been  described  with  an  emphasis 
on  phenomenological  language. 

A  major  objective  of  statistical  mechanics  is  to  predict  the  properties  for 
real  fluids  at  high  densities  from  the  known  properties  of  the  dilute  gas.  To 
accomplish  this,  it  is  usual  practice  to  reduce  exact  theories  (such  as  that 
associated  with  Eq.   (2))  to  theories  in  which  interactions  among  several 
particles  are  pairwise  additive.     The  properties  of  the  dense  fluid  can  then  be 
formulated  in  terms  of  the  detailed  structure  of  the  potential  functions  which 
describe  the  forces  between  pairs  of  particles.     Considerable  progress  has  been 
realized  using  this  approach.     Relevant  to  our  work  are  the  expansions  in 

density  based  on  Percus-Yevick  (PY)  theory  and  expansions  in  reciprocal  temper- 

19  20 
ature  using  the  Zwanzig  theory      as  modified  by  Barker  and  Henderson       (ZBH) . 

These  have  been  tried  for  several  potential  models  including  the  hard  and  soft 

spheres,  the  square  well,  and  the  Lennard-Jones  potential.     The  ZBH  temperature 

perturbation  theory  appears  to  be  the  more  successful  when  compared  to  results 

of  computer  experiments,  particularly  at  liquid  temperatures  but  also  for  gases 

at  high  densities  and  at  high  temperatures. 


Though  the  ZBH  theory  is  a  physically  satisfying  approach  and  does  compare 

well  (except  at  low  liquid  temperatures)  with  results  of  "computer  experiments", 

we  note  several  practical  limitations.     The  most  serious  of  these  is  associated 

with  the  application  of  the  theory  to  real  fluids.     In  the  ZBH  theory  the 

perturbation  terms  are  obtained  as  an  expansion  about  the  hard  sphere.  Barker 

20 

and  Henderson  have  shown      that  the  theory  can  yield  quantitative  results  when 

the  sphere  diameter  is  expressed  in  terms  of  a  kind  of  Boltzmann  average  of  the 

molecular  separation,  where  the  average  is  taken  over  the  positive  energy  region 

of  the  pair  potential.     Sphere  sizes  calculated  in  this  way  are  also  contained 

in  Figure  3.     It  is  a  fact,  however,  that  the  shape  of  the  pair  potential  even 

21 

for  simple  systems  is  quite  uncertain.     In  this  connection  it  has  been  proven 
that  for  realistic  (non-monotone)  potentials  it  is  not  possible  to  obtain  an 
unambiguous  representation  of  the  pair  potential  from  second  virial  coef f icients,, 
Also  at  high  densities  the  sphere  volume  tends  to  affect  the  equation  of  state 
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properties  somewhat  like  an  excluded  volume,  so  that  ambiguities  in  the  sphere 
diameter  are  amplified  in  their  effect  on  the  equation  of  state  at  high 
densities.    We  shall  present  an  example  of  this  below.     Thus,  though  the  ZBH 
theory  is  theoretically  satisfying,  its  application  would  seem  to  be  limited  to 
situations  where  the  pair  potential  is  known,  as  is  the  case  for  "computer 
experiments".     The  fact  that  the  use  of  this  equation  of  state  to  produce  thermo- 
dynamic tables  requires  a  complicated  numerical  integration  poses  a  second 
limitation  to  the  ZBH  approach.     Because  of  this,  results  have  so  far  been 
obtained  as  a  theoretical  end  in  themselves  and  are  therefore  of  somewhat 
limited  utility  to  the  engineer  or  scientist  who  desires  a  simple  analytic 
representation  of  the  equation  of  state  as  a  predictive  tool. 

It  should  be  appreciated  that  the  hard-sphere  diameters  are  only  conve- 
nient artifices.     Physical  interpretation  is  meaningful  only  in  the  context  of 
the  particular  overall  theory.     However,  since  the  ZBH  and  the  present  theory 
are  quite  sensitive  to  the  hard-sphere  diameter,  agreement  between  them  at 
least  to  within  several  percent  for  a(T)  would  be  necessary  for  the  two 
approaches  to  yield  comparable  equation  of  state  properties  at  high  densities. 

3.       Comparisons  of  the  HS  equation  of  state  with  experiment  and  with  other 
theories . 

In  this  section  the  HS  equation  of  state  is  used  to  calculate  compress- 
ibility factors  for  the  real  fluids  argon  and  molecular  nitrogen.     The  results 
are  compared  with  PVT  experimental  data,  and,  in  the  case  of  argon,  with 
results  of  "computer  experiments"  and  with  ZBH  theory.     The  experimental  second 
virial  coefficients  used  as  input  data  to  obtain  the  required  parameters  are 
smoothed  values  calculated  from  model  pair  potentials  as  determined  from  the 
experimental  second  virials. 

The  Figs.   4-8  contain  isotherms  calculated  for  argon  using  the  HS 
equation  of  state  and  the  ZBH  theory.     The  compressibility  factor  is  plotted 
versus  the  reduced  density,  p     =  Na3p,  for  isotherms  from  119. 8K  (=.8  critical) 
to  673. 15K.     Curves  labeled  #1  represent  the  compressibility  factors  predicted 
by  the  HS  equation  of  state.     The  curves  #2,  #3  and  #4  represent  results  of 
second  order  temperature  perturbation  ZBH  theory,   for  a  Lennard-Jones  gas  with 
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parameters  e/k  =  119. 8K  and  o  =  3.405  A  as  reported  in  references  22,  23  and  24, 
respectively:     the  curves  #2  and  #4  are  obtained  by  numerical  integration  of 
approximate  expressions  for  the  perturbation  terms  and  are  essentially  equiva- 
lent treatments;  the  curves  #3  refer  to  an  "exact"  calculation  of  the  second 
order  perturbation  terms  via  Monte  Carlo  techniques.     The  curves  #5  represent 
experimental  compressibility  factors  measured  for  argon.     The  data  points  on 
Figs.  4,  5  and  6  are  results  of  "computer  experiments"  for  the  Lennard- Jones 
gas  with  the  above  parameters,   the  circled  points  referring  to  Monte  Carlo 
results  and  the  boxed  points  to  molecular  dynamics  results.     The  "computer 
experiments"  could  be  uncertain  by  5  to  10%. 

Since  use  of  the  Lennard-Jones  potential  with  the  above  parameters  does  not 

25 

produce  second  virial  coefficients  which  fit  the  data  for  argon  below  200K  , 
Figs.  4-6  are  interesting  primarily  as  comparisons  of  theory  with  the  results  of 
the  "computer  experiments."     These  figures  show  that  the  present  theory,  as  well 
as  the  ZBH  theory,   are  only  qualitative  at  these  temperatures.     However,  the 
ZBH  curves  #3  yield  a  somewhat  closer  approximation  to  the  results  of  "computer 
experiments"  for  liquid  densities.     As  previously  stated,   the  Haar-Shenker 
equation  of  state  tends  to  degrade  at  low  temperatures.     But  it  is  apparent  from 
Fig.     4  that  the  theory  is  still  at  least  qualitatively  good  at  temperatures 
even  as  low  as  .8  of  the  critical  temperature.     In  fact  it  is  only  for  curves  #3 
(which  involve  extensive  numerical  calculations  to  evaluate  the  second  order 
perturbation  terms)  that  the  temperature  perturbation  results  are  a  significant 
improvement  over  the  HS  equation  of  state. 

An  explanation  for  the  fact  that  the  "wrong"  potential  gives  the  correct 

results  in  the  ZBH  theory  for  argon  has  been  offered  by  Barker,  Henderson  and 
26 

Smith.        They  argue,  that  the   (12,6)  pair-potential  with  the  above  parameters 
happens  to  be  an  "effective  potential"  that,   to  first  order  for  argon,  accounts 
for  high  density  non-additive  effects. 

The  Figs.   6-8  include  the  temperature  region  for  which  the  pair  potential 
used  produces  a  good  fit  to  the  experimental  second  virial  data.  Comparisons 
of  our  results  with  experimental  PVT  measurements  for  these  temperatures  are 
therefore  more  meaningful.     The  HS  equation  of  state   (curves  #1)   tends  to  follow 
the  PVT  experimental  data  (curves  #5)   fairly  closely  at  the  lower  densities,  up 
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to  about  p    =  0.6  (which  is  approximately  twice  the  critical  density).  At 
higher  densities  it  tends  to  yield  values  that  are  low,  but  in  most  cases  only 
by  less  than  5%.     The  results  are  roughly  comparable  to  those  for  the  more 
complicated  ZBH  theory.     The  ZBH  theory  (curves  #2,  //3,  #4,),  however,  tends  to 
underestimate  at  densities  near  the  critical  and  to  overestimate  at  high 
densities.     The  results  of  the  "computer  experiments"  seem  to  scatter  among  the 
various  theories.     Thus,  in  Fig.  6.    at  p*  =  .75,  the  Monte  Carlo  results  tend 

■k 

to  favor  the  HS  theory,  at  p    =  .9,  they  favor  the  ZBH;  while  at  p*  =  0.55  they 
fall  between  the  results  of  the  two  theories,  in  fact,  almost  on  the  PVT 
experimental  curve  #5. 

We  have  stated  above  that  uncertainties  in  the  pair  potential  for  real 

fluids  limit  the  utility  of  the  ZBH  theory  as  a  tool  for  predicting  the  equation 

of  state  properties  for  such  fluids.     To  illustrate  this  we  compare  the  equation 

of  state  of  argon  for  a  particular  isotherm  with  the  ZBH  and  HS  theories,  in 

which  different  inverse  power  representations  for  the  pair  potential  are  used, 
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each  of  which  produces  an  equivalent  fit  to  the  experimental  second  virials.  In 

Figs.     10  and  11  equation  of  state  results  are  presented  for  two  representations 

25 

for  the  pair  potential  for  argon:     the  (18,6)  with  the  parameters 

e/k  =  160.87  K 
a  -  3.261  A  . 

as  curves  #1;  the  (12,6)  with  parameters  given  earlier,  curves  #2.     Fig.  10 

includes  the  results  for  the  ZBH  theory;  Fig. 11  those  of  the  HS  equation.  Both 

figures  refer  to  the  isotherm  239. 8K.     In  these  figures  the  compressibility 

factor  is  plotted  against  the  density  in  amagats.     All  the  numerical  data  for 
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Fig. 10  were  furnished  by  Toxvaerd.         The  results  for  Fig. 10  show  the  two  curves 
near  coincidence  up  to  a  density  about  2/3  the  critical  (critical  density  for 
argon  =  300  amagats)  but  sharply  divergent  at  higher  densities.     By  comparison, 
the  results  for  the  HS  theory  in  Fig. 11  are  relatively  potential  independent. 
The  reason  for  this  independence  is  obvious.     The  HS  equation  depends  only  on 
the  experimental  second  virial  coefficient  and  its  first  derivative.     Since  the 
same  values  for  the  second  virial  coefficient  are  obtained  for  either  potential 
for  a  range  of  temperatures  near  this  one  (T        =  2.0),  the  calculated  second 
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virials  and  their  first  derivatives  are  equivalent.  Since  the  details  of  the 
repulsive  branches  of  the  potentials  are  different,  however,  the  calculation  of 
the  hard  sphere  diameter  according  to  the  prescription  of  Barker  and  Henderson 
yields  different  results  for  the  two  potentials.     This  difference,  in  turn, 
produces  a  difference  in  the  equation  of  state  predicted. 

Finally,  in  Figs.  12-15  we  present  results  of  using  the  HS  theory  for  calculat- 
ing the  equation  of  state  of  molecular  nitrogen.     The  compressibility  factor  is 
plotted  versus  the  fluid  density  in  amagat  units.     Here,  the  experimental 
second  virial  coefficients  were  represented  by  smoothed  tables  calculated  using 
the  Lennard-Jones  pair  potential  with  parameters  e/k  =  95.781K,  and  o  =  3.712  A. 
The  curves  #1  refer  to  the  present  theory;  the  curves  #2  to  experimental  PVT 
measurements.  The  latter  extend  to  10,000  atm  in  Fig.  15.     The  sphere  diameters 
are  obtained  from  Fig.   3  with  m=10.     The  comparison  of  theory  with  experiment 
for  nitrogen  is  quite  similar  to  that  for  argon.     As  with  argon,   the  HS  theory 
is  in  good  agreement  with  experiment  at  low  densities,  but  at  the  higher 
densities  the  experimental  PVT  isotherms  tend  to  be  slightly  steeper. 

4.  Summary. 

Based  on  the  virial  expansion  and  on  the  behavior  of  the  virial  coefficients 
at  high  temperature,  Haar  and  Shenker  derived  a  quantitative  yet  simple  equation 
of  state  which  is  valid  for  real  fluids  over  a  density  range  from  the  dilute 
gas  to  densities  approaching  that  of  the  solid  at  temperatures  above  twice 
critical,  and  which  requires  only  a  knowledge  of  the  second  virial  coefficient 
and  its  first  derivative  at  each  temperature.     This  equation  of  state  is  much 
simpler  than  the  ZBH  temperature  perturbation  theory  and  furthermore  does  not 
require  reference  to  the  precise  details  of  the  pair  potentials.     The  equation 
proved  to  be  quite  successful  in  comparisons  with  experimental  data.  These 
comparisons  were  naturally  carried  out  for  ordinary  temperatures  since  experi- 
mental data  exist  only  for  such  temperatures.     Figures  2  and  3  contain  such 
comparisons  for  argon  and  nitrogen. 

Through  the  use  of  the  intermolecular  potential  function,  the  HS  equation 
of  state  can  be  used,   in  a  very  simple  and  straightforward  manner,   to  extra- 
polate PVT  data  in  both  the  temperature  and  density  directions.     Thus,  the  inter- 
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molecular  potential  is  used  to  extrapolate  the  second  virial  coefficient  as  a 
function  of  temperature  and  this,  in  turn,  is  used  with  the  HS  equation  to  cover 
all  fluid  densities.     This  approach  is  therefore  particularly  well  suited  for 
aerodynamic  calculations  since,  as  already  mentioned,  these  require  temperature 
extrapolations  of  up  to  a  factor  of  ten  for  densities  up  to  those  which,  at  low 
temperatures,  correspond  to  the  liquid.     Because  of  the  decrease  in  the  attrac- 
tive contribution  to  the  second  virial  coefficient  with  temperature,  this 
equation  of  state  should  improve  with  increasing  temperature.    As  a  result,  the 
comparisons  made  at  ordinary  temperatures  should  easily  be  sufficient  for 
estimating  the  expected  adequacy  of  the  theory  at  high  temperatures.   Since  this 
equation  of  state  depends  only  on  the  second  virial  coefficient  which,  in  turn, 
is  determined  once  the  intermolecular  potential  is  known,  it  becomes  possible 
to  develop  an  entire  PVT  surface  given  this  intermolecular  potential  function, 
or  equivalently ,  given  a  sufficient  (generally  small)  amount  of  low  density  PVT 
data  for  the  substance  at  ordinary  temperatures,  from  which  data  the  inter- 
molecular potential  can  be  obtained. 

The  derivation  of  the  equation  of  state  is  based  on  an  expansion  in 
density,  where  the  reference  state  is  a  gas  of  hard  spheres.    We  have  presented 
a  plausibility  argument  which  indicates  that,  when  the  sphere  diameter  is  chosen 
appropriately,  the  terms  that  account  for  the  differences  between  the  properties 
of  the  actual  fluid  and  those  calculated  for  a  fluid  of  hard  spheres  are  sharply 
attenuated  at  temperatures  above  the  liquid-vapor  critical  temperature,  for 
densities  up  to  that  of  the  solid. 

It  has  been  shown  that  the  temperature  dependence  of  the  equation  of  state 
at  temperatures  above  the  critical  is  determined  by  two  parameters  which  depend 
on  temperature  and  that  these  parameters  can  be  obtained  from  the  second  virial 
coefficient  and  its  first  temperature  derivative. 

III.     THE  EXTENSION  OF  THE  HAAR-SHENKER  EQUATION  OF  STATE  TO  MIXTURES 

Haar  and  Shenker  developed  their  equation  of  state  for  use  with  pure  sub- 
stances.    Our  needs  are,  of  course,  for  a  theory  applicable  to  mixtures  since 
at  aerodynamic  temperatures  even  pure  nitrogen  becomes  a  mixture  as  a  result  of 
dissociation  and  ionization.     There  are  several  ways  in  which  a  theory  for  pure 
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fluids  can  be  extended  to  mixtures.     Perhaps  the  simplest  are  what  have  been 
called  the  one  and  two  fluid  van  der  Waals'  theories.     In  the  one  fluid  model 
the  mixture  behaves  as  if  it  were  a  single  component  fluid  (which  can  be  called 
the  equivalent  fluid)  with  any  parameters  used  to  describe  the  fluid  being 
averaged  over  the  composition.     This  averaging  is  carried  out  by  taking  the 
parameters  associated  with  the  individual  constituents  and  suitably  weighting 
them  to  the  extent  that  the  constituents  are  present  in  the  fluid.     For  example, 
a  fluid  made  up  of  a  mixture  of  hard  spheres  would  be  described  as  a  one  fluid 
van  der  Waals'  model  by  the  equations  associated  with  a  single  component  hard 
sphere  fluid  but  with  the  single  relevant  parameter  (that  associated  with  the 
molecular  diameter)  averaged  over  the  composition.     A  natural  way  of  doing  this 
is  to  take  for  the  hard  sphere  volume 


where  a      is  the  diameter  of  the  molecules  of  the  equivalent  fluid  and  a.. 

eq  th  ij 

the  diameter  for  the  interaction  between  a  molecule  of  the  i      species  and  one 

of  the         species  in  the  actual  fluid.     For  hard  spheres  a..  —  =  (a.  +  a.). 

ij      2      1  3 

The  equation  of  state  for  this  mixture  is  easily  derived.     From  the  virial 

theorem,  it  is  possible  to  derive  a  general  equation  of  state  for  a  fluid  in 

28 

terms  of  the  distribution  of  pairs  of  particles  in  the  fluid.        Because  of  the 
abruptness  of  the  hard  sphere  interaction,  this  equation,  for  a  one  component 
hard  sphere  fluid,  reduces  to 


where  g(j)  is  the  probability  that  a  pair  of  molecules  will  be  found  a  distance 
a  apart.     For  a  mixture  of  such  spheres,  the  equation  of  state  becomes 


o3      =  T      X,  X 
eq  i 


pa3  g(a) 


(10) 


X.  X.  a3, 
i    1  ij 


(11) 
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ttl 

where  g. .    (a..)  is  the  probability  that  a  pair  of  molecules,  one  of  the  i 
1J       1J  th 

species  and  one  of  the  j  species,  are  separated  by  a  distance  a .  In  the  one 
fluid  model,   then,   (11)  is  replaced  by  (10)  in  which 

§ij   (o..)  =  g  (a) 
and  (12) 

o3  =  y  x.  x.  o?. 


This  model  is  consistent  with  the  results  of  the  density  expansion  in  the 
Ursell-Mayer  virial  equation  of  state. 

A  number  of  models  other  than  (12)  can  be  devised  for  using  (10)  in  place 
of  (11)   for  a  mixture.     Each  of  these,   though  reasonable,  does  not  lead  from 
(11)   to  (10)  in  a  natural  way.     One  might,   for  instance,  average  the  diameter 
rather  than  the  volume  so  that 

a  =  y   X.  X.  a. .  =  y  X.  a. .  since  a..  =  -  (a  .  +  a  ) 
^     i    j     xj  x    xi  ij      2  v  xx  23 

i.  j  i 


One  might  also  average  the  volumes  over  like  species  only,  i.e.   a3  =  E  X_^  CT|^- 
We  shall  consider  (12)  as  the  only  reasonable  model,  particularly  since  it  is 
the  only  one  that  leads  to  results  which  are  satisfactory. 

A  two  fluid  theory  can  be  obtained  from  the  approximation 


,     (0     )  =  i U  .   (a.  .)  +  g.  .   (a.,)  I 

'i:    iy     2L  11     11       j  J  J 


If,  at  the  same  time,  one  uses  the  fact  that,   for  hard  spheres  a . .  =  r  (o.+  a.) 

ij       2      x  j 

(11)  becomes,  for  this  model, 
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RT  ~  1  +  ~  P    L  Xi  °-  8±i  (o.  . )  (13) 
i 

In  order  for  this  to  be  of  the  pure  fluid  form,  it  must  be  assumed  that  g  is 
the  distribution  function  for  a  pure  fluid  of  hard  spheres  of  volume  a.  given 
by 

2 


Thus,  if  this  were  a  binary  system,  it  would  appear  to  be  made  up  of  two  pure 
fluids  of  different  diameters,  namely  a3  =        (Tjj  +  X2  a32  and  a3  =  X1  ar^  + 
X2  °22'     Hence  the  name  two  fluid  theory. 

A  fluid  of  hard  spheres  is  a  highly  idealized  model  for  an  actual  fluid  so 

might  be  thought  to  be  quite  useless  for  testing  theories  of  fluids.     There  are, 

however,  properties  for  fluids  of  hard  spheres  as  calculated  by  computer 

simulation  methods.     If  such  results  obtained  for  hard  sphere  mixtures  are  taken 

as  "experimental"  data,  it  is  then  reasonable  to  compare  them  with  (12)  and  (13) 

to  see  which  is  the  better  approximation  to  a  mixture  of  hard  spheres.  This 

29 

has  been  done  by  Henderson  and  Leonard      who  found  the  one  fluid  theory  to  be, 
by  far,  the  superior. 

In  a  later  paper,  these  same  authors  carried  out  a  similar  comparison"^  for 

a  fluid  whose  molecules  interact  in  accordance  with  a  Lennard- Jones  (12,6) 

potential  function.     Such  a  potential  contains  both  attraction  and  repulsion, 

and  is  often  characterized  by  a  molecular  diameter  and  by  a  potential  well  depth. 

The  one  fluid  model  now  follows  from  ea3  =    E    X.  X.  e..  a3...     Comparisons  were 

±a     i    J     ij  ij 

made  for  an  equimolar  mixture  and  the  one  fluid  model  was  found  to  be  superior, 
particularly  in  the  prediction  of  the  excess  free  energy  and  heat  of  mixing. 
The  intercomparison  for  the  excess  volume  of  mixing  was  somewhat  ambiguous, 
however. 
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It  would  thus  seem  clear  that  the  most  reasonable  way  in  which  a  theory  for 
a  pure  fluid  can  be  applied  to  that  for  a  mixture  is  via  a  one  fluid  model. 
This  we  have  chosen  to  do  for  the  HS  equation  of  state.     This  simplifies  our 
task  since  we  need  only  write  down  the  pure  fluid  form  of  the  equation  of  state 
which  we  have  already  derived  and  interpret  the  parameters  in  the  form  of  the 
one  fluid  theory. 

By  (5),  the  HS  equation  of  state  for  a  pure  fluid  is  written 


pv  =  i±v+£  +  4/5  _  A  (14) 


y 

(l-y) " 


RT    „     ,3  \b 


where  y  =  .  Then  Eq .  (14)  becomes  an  equation  of  state  for  a  mixture 
according  to  the  one  fluid  model,  if 

b  =    Y\     X.  X.  b. .  (15) 
i  J      1    2  13 


b  =  t  x.  x. 

4      A         1  J 


B.  (16) 


It  should  be  remembered  that  Eq.  (9)  must  be  solved  for  b„  for  each  of  the 
interactions  included. 


IV.     THERMODYNAMIC  FUNCTIONS  IN  THE  ONE  FLUID  HS  MODEL 


In  this  section  of  the  report  we  shall  be  concerned  only  with  the  real  gas 
contributions  to  the  various  thermodynamic  properties.     The  ideal  gas  contri- 
butions and  certain  precautions  required  in  their  calculation  are  contained  in 
reference  (3) . 

A.       HELMHOLTZ  FREE  ENERGY,  CHEMICAL  POTENTIAL,  AND  ACTIVITY  COEFFICIENTS. 

The  combination  of  the  equation  of  state  (14)  and  the  one  fluid  model,  (15) 
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and  (16),  is  sufficient  for  the  calculation  of  the  equation  of  state  of  a  real 
gas  mixture  at  high  densities.  With  such  an  equation  of  state  it  then  becomes 
possible  to  calculate  the  real  gas  part  of  the  Helmholtz  free  energy  from 


from  which  the  chemical  potential  can  be  calculated  using 


ui      \dnj  T,  p,  n_. 


when  n.  refers  to  all  n.,  j  ^  i. 
J  J 


When  the  ideal  and  real  parts  of  the  chemical  potential  are  separated  it 
becomes  possible  to  identify  the  activity  coefficient  for  each  species.  By 
properly  combining  the  activity  coefficients  for  each  species  taking  part  in  a 
chemical  reaction,  one  can  define  an  effective  activity  coefficient  for  that 
chemical  reaction.     This  effective  activity  coefficient  directly  modifies  the 
equilibrium  constant  to  produce  the  effect  of  non-ideality  on  the  chemical 
reaction.     To  see  how  this  goes,  let 


A  „. 


where  is  the  ideal  gas  part  of  the  chemical  potential  and  Ay^  the  real 

part.  It  follows,  then,  that  the  activity  coefficient  for  this  species  is 
given  by 


Ay. 

w  =  lny± 


The  effective  activity  coefficient  for  species  i  in  a  chemical  reaction  with 
stoichiometric  coefficients  v  . ,  is  then  given  by 
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which  is  equivalent  to 


y]  =  exp  [(E  V. .  Au.  -  Ay  ) /RT] 
j  3 


(17) 


We  shall  now  proceed  to  the  detailed  calculation  of  the  free  energy  and 
chemical  potential  for  the  HS  equation  of  state. 

If  at  each  volume  the  reference  state  for  the  free  energy  is  taken  as  that 
of  the  ideal  gas  at  that  volume  then 

V 

A(V)  =  A(«)  -f  PdV 


A(0)(V)     =  A(o)(~)  -  /  P 


<•>  dV 


But,  since     A  (oo)        =  (oo)  it  follows  that 


A(V)  -  A(o)(V)  =  -   /     (P  -  P(0)) 


dV 


nZRT      j  ^(o)      nRT      ,         „  .     ^,  r    ^  , 

Since  P  =  — —  and  P        =  — — ,  where  Z  is  the  compressibility  factor  and  n  the 


total  number  of  moles,  it  follows  that 

A(V)  -  A(o)(V)  =  -  nRT  f  {1-1)  ^ 

OO 

The  real  part  of  the  chemical  potential  is  then  given  by 
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RT 


y  .- 
i 


(o) 
i 


RT 


The  evaluation  of  the  common,  integral  is  obviously  necessary  to  progress  in  the 
derivation  of  both  the  free  energy  and  the  chemical  potential. 

According  to  the  HS  theory, 


Cl-y)3  b 


>^2-  2>T+4)       +  (|  -  1)  4y 

Cl-y)3  b 


Since  y  =  7-=:,  —  =  -  &  .  Thus 
■»       4v'     V  y 


Cl-y)3  „        b  1  ' 


In  1^,  let  1-y  =  x  so  that  dy  =  -dx  and 


h  .j-^'a-x)2  -2(i-x)+4  dx 


1  x- 


1  x 


Since,  obviously, 
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I2  =  -  4y(E  -1) 


then 


/ 


CZ-l)  f-  =  -  |  [Ojry)2  -1]  +  ta  (l-y)  -  4y  (5  -1)  Q8) 


so  that 


A(V)  -  A(Q)(V) 
RT 


%  [(IT>2  -U  -n  to  (l-y)  +  4ny(?  -  1) 


The  non-ideal  contribution  to  the  chemical  potential  is  then  given  by 


Ay 
RT 


The  first  term  has  already  been  evaluated.  The  second  is  obtained  by  differ- 
entiating (18).  Thus 


3  (-2)(-l)       1  B      .J   9y       .     3  .B, 

2     "     ,3    -lTy  -  X)     an"."  ^  to.^ 

(l-y)  7  )       l  i 


At  this  point  use  must  be  made  of  the  one  fluid  model.  Thus 


1       3b  1 


3n.       4V     3n.  4V 
l  l 


2  E  n.  b.. 
(E  n.)2 


2  En.  n.  b.,  . 


(E  n.) 


1_  2 
4V  n 


En.  b.  . 


-  b 
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Also 


9  1     3B        B_  3b 

"ST  CB/b)  "  b     3n±      b2  3n± 


3B  9b 

Since,  obviously,  - —    will  have  the  same  form  as  —    but  with  B  replacing  b, 

dn^  dn^ 

we  can  immediately  write  down 


dn.   <V  2 


n  b 


I  n.  B. .  -  r  E  n  b 
I    ]    lJ      b  .    j  ij 


Combining  all  results  obtained  to  here  yields 


^i  _  3 
RT  2 


£n  Ci-y)  +4  C-  -  Dy 

b 


3      1-y  b 


d-y) 


22 
b 


I  n  .  b  .  . 
....I  jJ 


,  2 
+  4y  5 


"En.  B. . 


*2  En.  b.  . 
B  j  ij 


For  simplicity  we  write  this  as 


En.  b.  . 

Ay.  j  u 

RT"      =  yo  +  yl   nb  


E  n.  B.  . 
+  y2  ^-^b— 


where 
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yo  =  2 


-1 


Xi-y)' 


-  In  (1-y)  - 


^i-  2?\—  3  +  1=7  "  4! 

1  [Cl-y)J         7  j 


P2  =  8y 


Substitution  in  (17)  then  yields 


=  exp 


-    cj    u    +  — —      E  n./Ev..b,n-b.„) 
i    o        nb      %      %\.     ij     j£  i^/ 


+  ^  n„  / E  v..   (B.„-  B.n)\  /RT 

nb         t    I  L  xj 


C19) 


where  -oj.  =  E  v..  -1  is  the  net  decrease  in  the  number  of  particles  in  the  i 
reaction  defined  earlier.    With  this  expression  for  y  '  it  becomes  possible  to 
calculate  the  effect  of  the  HS  theory  on  the  equilibrium  composition  of  the 
mixture . 

B.      ENTROPY,  ENTHALPY  AND  GIBBS  FREE  ENERGY 

It  is  now  possible  to  calculate  all  of  the  thermodynamic  properties  predict- 
ed for  the  HS  model.     The  entropy  follows  from 


so  that  the  real  gas  part  is  given  by 
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AS 
R 


VIA 
3T 


AA  T  _3  /1A 
RT  ST  \RT 


But 


3^    +4  M 


+  4ny  ty^). 


Now 


and  3T  Cb}  "    b     3T        bZ  BT 


so  that  -r=      TTT  = 


1A 
T     \  RT 


l(l-y) 


+  J-  -4  1  »  |B 
3        1-y        |     b     ST        b  3T 


Thus 


_S_ 

R 


3n 
2 


c-U2  -1 

1-y 


-r.  in 


CL-y)  +  4  n  y  G  -D 


(Cl-y)3      ^       i    b  3T 


4yT  _3B 
b  3T 


3b  3B 

Trrr  and  -r=  are  expected  to  be  quite  small  because  of  the  high  temperatures  of 

dl  dl 

interest.  It  is  therefore  possible  that  the  last  two  terms  might  be  negligible 
with  respect  to  the  first  three.     This  would  need  to  be  examined,  however. 


This  expression  for  the  entropy  and  the  earlier  one  for  the  equation  of 
state  can  be  used  to  calculate  the  enthalpy  and  Gibbs  free  energy.  Specific 
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heats  can  be  calculated  from  these  functions  either  by  numerical  differentiation 

or  by  means  of  algebraic  expressions  which  can  be  derived  by  differentiation  of 

3 

the  above  expressions  . 

V.     NUMERICAL  METHODS 

The  primary  step  in  calculating  the  thermodynamic  properties  consists  of 
the  determination  of  the  equilibrium  composition  of  the  mixture.     This  requires 
the  solution  of  equations   (1)  subject  to  the  constraint  that  the  total  mass  of 
each  nuclear  type  be  conserved.     The  effective  activity  coefficients,        ,  are 
given  by  equation  (19)  for  the  contribution  of  the  HS  model  to  the  equilibrium 
constant  and  in  appendix  B  of  reference  (3)  for  that  of  the  Debye-Huckel  theory 
of  ionic  solutions.     There  are  two  lower  level  approximations  to  the  HS  model 
which  we  have  already  considered  in  earlier  work.    We  have  already  used  the 
ideal  gas  approximation,  in  which  all  y  '  are  taken  equal  to  unity.     We  have 
also  made  use  of  the  second  virial  coefficient  approximation  for  which  the  y^* 
of  the  HS  model  are  replaced  by  the  term  linear  in  the  density  in  their  density 
expansions  (which  term  appears  in  appendix  B  of  reference  (3)). 

The  method  used  for  the  solution  of  the  equations  for  the  concentrations  in 

both  the  ideal  gas  and  second  virial  gas  has  been  described  in  some  detail  else- 

3 

where.  The  procedure  used  here  was  required  to  be  applicable  to  all  three 

approximations,   i.e.   the  ideal  gas,  second  virial  and  HS.     Although  based  on 
these  earlier  methods , there  was  a  considerable  modification  of  certain  details 
of  the  original  approach  to  allow  for  greater  flexibility  in  obtaining  solutions. 
In  order  to  permit  reference  species  to  be  chosen  at  will  simply  through 
modification  of  the  input  data  (as  described  in  reference  (3)),  the  calculation 
of  the  equilibrium  constants  was  made  part  of  the  computer  program.     This  led  to 
the  discovery  of  an  error  in  the  equilibrium  constant  used  for  0^  in  the 
previous  NBS  tables.     Our  correction  of  this  error  caused  the  search  method  used 
for  finding  the  electron  concentration  to  be  unstable  at  low  temperatures.  A 
not  insignificant  amount  of  time  was  spend  in  isolating  this  problem  and  in 
correcting  it.     The  problem  was  associated  with  orders  of  magnitude  increase  in 
the  concentration  of  0^  which  occurred  after  the  correction  was  inserted.    As  a 
result,  the  electron  concentration  now  became  orders  of  magnitude  smaller  than 
the  concentration  of  0_  and,  in  fact,  was  now  calculated  as  a  small  difference 
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between  the  concentrations  of  NO  and  CL.  This  often  led  to  negative  interim 
values  for  the  electron  concentration  during  the  non-linear  search  procedure. 
Since  the  electron  concentration  has  only  a  minor  effect  on  the  values  obtained 
for  the  other  reference  species,  the  problem  was  solved  by  instituting  a  grid 
search  for  the  electron  concentrations,  over  positive  values  only,  whenever  a 
negative  guess  value  was  obtained  for  that  concentration,  leaving  all  other 
concentrations  fixed  at  their  current  values. 

At  low  temperatures,  it  was  also  necessary  to  produce  a  very  strong 
"damping"  of  the  electron  concentration  by  means  of  a  q  value      smaller  than 
(but  in  the  neighborhood  of)  unity.     This  reduced  the  "natural"  excursions  in 
the  electron  concentration  from  iteration  to  iteration. 

Initially,  precisely  the  same  search  procedure  was  used  for  the  HS  model 
as  was  used  earlier  for  the  ideal  and  second  virial  gases  except  that  the 
for  the  HS  model  now  had  to  be  computed  at  each  iteration.     Problems  arose, 
however,  because  these  y.'  became  quite  large  for  certain  species  at  the  highest 
densities.     In  fact,  y.'  values  approaching  10    were  encountered.     This  produced 
instabilities  in  the  search  procedure  used.     This  problem  was  solved  by  what 
might  be  called  a  dual  level  search  procedure.     In  this  procedure  an  initial  set 
of  y  1  were  computed  based  on  the  initial  guesses.     These  y.1  were  maintained 
constant  until  a  set  of  concentrations  was  obtained  which  satisfied  the  mass 
balance  equations  for  the  gas  with  these  initial  y  '  values.    With  these 
concentrations,  a  new  set  of  y^'  was  computed.     These  new  y.'  were  now  held 
constant  and  a  second  set  of  concentrations  obtained  which  satisfied  the  mass 
balance  equations  for  this  second  set  of  y.'  values.     The  procedure  was  carried 
out  repeatedly  until  the  largest  change  in  the  y^'  when  recomputed  was  less  than 
a  given  tolerance.     To  reduce  the  possibility  of  producing  numerical  instabilit- 
ies, only  a  fraction  of  the  change  calculated  for  each  y.1  was  used  in  any  new 
iteration  and,  in  any  case,  the  y  '  were  not  allowed  to  change  by  more  than  some 
arbitrary  factor.     In  order  to  reduce  the  computer  time  required,  a  coarse 
tolerance  was  placed  on  the  values  of  the  concentrations  accepted  as  correct  for 
the  first  few  sets  of  values  (generally  incorrect)  obtained  for  the  y  1   in  this 
procedure. 

The  present  computer  program  has,  to  a  large  extent,  retained  the  ability 
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of  the  earlier  program  to  obtain  solutions  essentially  independent  of  the 
initial  guess  values  for  the  concentrations  of  the  reference  species.  The 
electron  concentration  is  perhaps  an  exception  to  this  but  then  mainly  only 
when  the  computation  on  an  isotherm  begins  at  quite  high  densities. 

Since  the  few  problems  which  we  have  encountered  with  obtaining  solutions 

must  depend  on  the  specific  search  procedure  used  by  us,  it  would  be  useful  to 

study  the  process  of  solution  for  other  search  procedures.     Because  of  its 

simplicity,  we  would  place  particular  emphasis  on  a  study  of  the  direct  search 

31 

method  of  Hooke  and  Jeeves 

Once  the  species  concentrations  are  obtained  in  the  above  manner  (with  the 
mass  balance  equations  satisfied  to  the  tolerances  specified  for  each  reference 
species) ,   the  thermodynamic  properties  of  the  mixture  can  be  calculated  in  a 
straightf orward  manner. 

VI .  RESULTS 

In  this  section  we  discuss  a  number  of  special  features  of  the  results 

obtained,  reserving  our  discussion  to  those  results  which  illustrate  the  effects 

of  density.     We  compare  results  obtained  by  us  for  the  ideal  gas,   second  virial 

coefficient  and  HS  models  with  each  other  and,  where  appropriate,  with  the 

results  of  the  GB  extrapolation.     We  shall  also  compare  the  concentrations 

predicted  for  various  species  among  the  three  approximations  used  by  us.  The 

8 

Grabau-Brahinsky  model  does  not  include  the  calculation  of  the  species  concentra- 
tions so  cannot  be  included  in  that  part  of  the  discussion. 


Certain  interesting  features  of  the  predicted  dependence  of  the  concentra- 
tions on  density  are  also  discussed.     Of  particular  interest  is  the  demonstration 

32 

of  the  "strength"  of  Le  Chatelier's  principle 

A.       THE  DEPENDENCE  OF  SPECIES  CONCENTRATIONS  ON  DENSITY. 


The  formalism  which  we  have  developed  has  the  capability  of  predicting 
density  effects  on  species  concentrations.     This  ability  is  important  for  two 
reasons.     It  is  obviously  important  when  the  concentrations  of  the  species  them- 
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selves  are  needed,  as  might  be  the  case  when  specific  species  have  particularly 
interesting  properties.     Examples  might  be  the  electrically  charged  species 
(which  obviously  affect  the  electrical  conductivity  very  strongly)  and  any 
species  which  radiate  in  a  useful  part  of  the  spectrum.     Being  able  to  predict 
density  effects  on  concentrations  is  also  important  in  the  calculation  of  the 
density  effects  on  the  thermodynamic  properties.     The  expressions  for  the  thermo- 
dynamic properties  of  the  reacting  mixture  are  sums  over  concentrations  times 
properties  of  individual  species  plus  sums  of  products  of  concentrations  times 
properties  of  pairs  of  species.     In  our  model,  the  latter  enter  into  the  calcul- 
ation of  the  average  (i.e.  one  fluid)  molecular  diameter.     The  thermodynamic 
properties  of  the  mixture  therefore  depend  on  temperature  and  pressure  through 
the  dependence  of  the  concentrations  on  these  state  parameters. 

Density  effects  on  the  thermodynamic  properties  of  air  have  generally  been 
calculated  by  others  on  the  assumption  that  the  species  concentrations  either 
do  not  change  with  density  or  simply  obey  the  ideal  gas  mass  action  law.  This 
was  essentially  the  assumption  of  Grabau  and  Brahinsky.     As  we  shall  show,  this 
kind  of  assumption  breaks  down,  particularly  at  the  highest  densities  where  the 
identities  of  the  major  species  in  the  mixture  change,  when  density  effects  are 
specifically  taken  into  account.     This  can  change  the  entire  character  of  the 
gas  being  studied. 

Care  must  be  taken  here  not  to  place  too  much  emphasis  on  the  actual 
numerical  values  obtained  by  us  for  the  concentrations.     The  species  concentra- 
tions at  high  densities  are  expected  to  be  sensitive  functions  of  the  assumptions 
made  in  any  calculation  of  this  kind.     In  our  model,  they  can  be  expected  to 
depend  very  strongly  on  the  intermolecular  forces  used  (especially  at  the  high- 
est densities)  as  well  as  on  the  method  of  calculating  the  hard  sphere  diameter 
at  temperatures  for  which  the  second  virial  coefficient  has  a  negative  slope. 
This  problem  is  associated  with  the  fact  that  the  effective  hard  core  potential 
produces  a  second  virial  coefficient  whose  slope  is  always  positive.  The 
relationship  of  the  actual  numerical  values  obtained  by  us  to  actual  air  might 
also  be  expected  to  depend  strongly  on  our  omission  of  such  species  as  ^0^,  C^t 
0^  and  3ome  of  the  molecular  ions. 

Because  of  Le  Chatelier's  principle,  errors  introduced  into  the  calculation 
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of  the  concentrations  for  any  one  of  these  reasons  should  be  less  than  that 
which  might  be  estimated  from  a  direct  calculation  of  the  apparent  effect  of  the 
omission.     This  principle  is  essentially  a  statement  of  the  competition  among 
the  various  chemical  reactions,  because  of  which  species  concentrations  "resist" 
change.     Thus,  suppose  a  correction  to  the  equilibrium  constant  in  a  certain 
reaction  produces  an  increase  in  the  concentration  of  a  particular  species. 
This  increase  will  cause  a  competing  reaction  containing  that  same  species  to 
move  in  such  a  direction  as  to  reduce  the  concentration  of  this  particular 
species.     The  overall  effect  of  the  original  correction  to  the  equilibrium 
constant  is  thereby  reduced.     An  example  of  this  is  described  below.  Because 
of  this  and  because  any  errors  resulting  from  our  various  assumptions  can  be 
expected  to  affect  the  different  species  perhaps  randomly  according  to  sign,  the 

ttl 

overall  effects  of  our  zero      order  assumptions  on  the  properties  computed  should 
actually  be  much  smaller  than  might  be  expected  from    an  examination  of  the 
effects  of  the  various  separate  approximating  assumptions.     This  error  reduction 
might  be  considered  to  be  a  decided  advantage  in  favor  of  the  use  of  a  micros- 
copic molecular  model  such  as  ours. 

The  operation  of  Le  Chatelier's  principle  can  be  seen  by  comparing  the 
change  in  various  concentrations  actually  obtained  in  the  calculation  when  the 
density  is  changed  at  constant  temperature  with  that  which  might  have  been 

-0) . 

expected  from  the  change  in  the  density  factor  (p/po)     1  in  equation  (1). 

Consider,  for  example,  going  from  log  p/pq  =  2.0  to  3.0  at  T  =  3000K.     For  the 

species  N„0,  a  concentration  enhancement  by  a  factor  of  1.6x10^  could  have  been 

4 

expected  based  on  the  value  y.     =  4.3x10  at  log  p/p    =  3.0  and  to.  =  0.5,  whereas 

X  ox 

the  actual  enhancement  obtained  involved  only  a  factor  of  1.7x10^.  There 
was,  therefore,  a  reduction  by  a  factor  of  100. 

Figures  16  and  17  contain  plots  of  concentrations  against  density  for  a 
number  of  species  for  the  temperatures  3000K  and  9000K.     The  concentrations  are 
those  which  were  calculated  using  the  full  density  effect  with  the  HS  equation, 
and  those  based  on  the  ideal  gas. 

The  dependence  on  density  of  the  concentration  of  oxygen  and  the  effect  of 
this  dependence  on  the  concentration  of  the  other  species  are  particularly 
interesting.     The  rapid  decrease  in  the  concentration  of  molecular  oxygen  at 
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the  highest  densities  is  very  dramatic  for  the  HS  gas,  especially  when  compared 
to  its  behavior  in  the  other  two  approximations.     As  a  result,  the  gas  at  3000K 
and  1000  times  normal  density  consists  mainly  of  ^  and  NO  rather  than  of  ^  and 
O2  as  predicted  in  the  ideal  gas  approximation.      Secondary  effects  from  this 
reduction  in  0^  concentration  can  be  seen  in  the  marked  reduction  in  the 
concentrations  of  0~,  0*  and  0    predicted  by  the  HS  theory  as  compared  to  the 
predictions  of  the  ideal  gas  approximation. 

The  reduction  in  0£  concentration  "results"  from  the  enhancement  in  the 
production  of  NO,  NO£  and  ^0  and  the  associated  requirement  for  the  production 
of  atomic  oxygen  for  the  formation  of  these  species.     The  cause  of  the  enhance- 
ment of  these  species  can  be  seen  in  Table  2  which  contains  the  y.'  values  for 
all  species  at  T  =  3000K  for  several  densities.     In  the  concentration  units  used 
by  us,  it  is  necessary  to  multiply  the  y      by  (p/p  )     1  to  obtain  the  full 
density  effect  along  an  isotherm.     As  is  pointed  out  below,  in  some  instances 
this  produces  enhancement  factors  of  10^  for  the  equilibrium  constant. 


An  interesting  effect,  which  might  be  called  a  second  order  effect  of  the 
reduction  in  0£  concentration,  was  seen  among  the  carbon  containing  compounds. 
According  to  the  y^'  values  of  Table  2}  the  concentration  of  C  should  drop 

drastically  with  increasing  density.     This  is  especially  true  when  the  factor 

-1  -3 

Cp/p  )      is  added  (bringing  in  an  additional  factor  of  10      at  the  highest 
o 

density).      Instead,  the  concentration  of  C  increases  slightly  with  density. 
This  comes  about  because  the  reduction  in  0^  concentration  causes  a  decrease  in 
the  production  of  CO2.     This  "frees"  carbon  atoms  which  become  available  for 
the  enhancement  of  the  CO  and  C  concentrations. 


Another  second  order  effect  is  that  of  the  slight  increase  in  electron 
concentration  with  increasing  density,  quite  the  opposite  of  the  ideal  gas 
behavior.     This  results  from  the  decrease  in  the  concentrations  of  O2  and  0 
which  results  in  additional  free  electrons.     The  increase  in  the  electron 
concentration,  in  turn,  results  in  a  drop  in  the  concentration  of  N0+,  the  main 
electron  producer,  causing  the  increase  in  electron  concentration  to  be  somewhat 
reduced  from  that  expected  purely  on  the  basis  of  reduction  in  O2,  another 
example  of  Le  Chatelier's  principle. 
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Those  species  for  which  y^'  has  the  value  of  unity  in  Table  1  are  either 

reference  species,  for  which  no  equilibrium  constant  is  needed,  or  are  ions  for 

which  no  estimates  of  virial  coefficients  were  available.     It  should  be  noted 

that  the  y^1  values  for  the  species  ^0  and  N0^  begin  to  approach  10^  at  the 

highest  densities.     This  effect  is  further  enhanced  by  the  factor  33    for  each 

-1/2 

arising  from  the  factor  (p/pQ)        .     The  overall  density  effect  for  these 
species  at  p/pq  =  1000    as  compared  to  p/PQ  =  1«0  is  then  such  as  to  multipy  the 
equilibrium  constant  by  a  factor  of  over  10^.     It  is  therefore  obvious  why  there 
is  such  a  strong  dependence  of  concentration  on  density. 

The  competition  at  the  highest  densities  among  NO,  NC^  and  ^0  (and 
especially  between  the  last  two)  is  particularly  interesting.     As  the  density 
increases,  all  three  species'  concentrations  increase.    At  3000K  this  occurs 
mainly  at  the  expense  of  molecular  oxygen.     As  the  density  increases,  however, 
the  enhancement  of  N^O  begins  to  proceed  at  such  a  pace  as  to  "require"  oxygen 
atoms  from  other  reactions  so  that,  ultimately,  the  increase  in  N^O  concentra- 
tion takes  place  at  the  expense  of  the  concentrations  of  NC^,  and  NO.  At 
9000K,  the  initial  enhancement  occurs  at  the  expense  of  the  oxygen  atom  con- 
centration but,  at  intermediate  densities,  produces  a  reduction  in  molecular 
oxygen  concentration.     Eventually,  the  increased  N£0  concentration  occurs  as  a 
result  of  a  reduction  in  NO  and  NO^  concentrations. 

The  absence  of  estimates  of  the  virials  for  interactions  involving  ionic 
species  obviously  leads  to  errors  at  the  highest  temperatures  where  charged 
species  become  non-negligible.     Although  these  charged  species  do  not  dominate 
in  our  approximation,  it  is  conceivable  that  they  might  become  major  constitu- 
ents through  enhancements  caused  by  large  y  1  values  which  might  be  obtained. 
This  effect  might  be  expected  to  be  smaller  for  the  ionic  species  than  it  was 
for  the  neutrals  NO,  N0„  etc.  for  two  reasons.  First  of  all,  even  at  high  temp- 
eratures, repulsion  between  ions  of  like  sign  is  reduced  by  the  effects  of 
attraction  between  those  of  opposite  sign.     Secondly,  our  inclusion  of  the  Debye- 
Huckel  limiting  law  (see  reference  (3))  already  Includes  part  of  the  interaction 
between  charged  particles. 
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B.       DENSITY  EFFECTS  ON  THE  EQUATION  OF  STATE. 

Table  3  contains  the  compressibility  factor  as  a  function  of  density  at 
selected  temperatures  for  all  three  approximations  as  well  as  for  that  of  Grabau 
and  Brahinsky.     The  ideal  gas  approximation  contains  no  correction  to  the 
equation  of  state  for  the  effect  of  density.     Since  there  is  a  "density  effect" 
in  the  law  of  mass  action  as  written  in  the  units  used  by  us  there  is  a 
variation  in  the  calculated  compressibility  factor  with  density.     By  our 
definition,  the  compressibility  factor  for  the  ideal  gas  is  simply  the  total 
number  of  moles  of  the  mixture.     According  to  the  equations  of  chemical  equilib- 
rium for  the  ideal  gas  as  written  in  our  units,  the  effective  equilibrium 

constant  for  a  particular  reaction  in  the  ideal  gas  approximation  is  given  by 

eff  .  nij  -co.. 

K.      =  K.   (T/T.)     1  (p/p  )     1  and  so  increases  with  density  for  those  reactions 

i  l         (J  o 

in  which  a  net  decrease  in  the  number  of  particles  results  in  the  production, 
from  the  reference  species,  of  the  species  associated  with  the  reaction.  Such 
a  reaction,  for  example,  is 

N20  =  N2  +  |  02 

where  two  molecules  are  produced  for  every  three  which  react.     On  the  other 
hand,  a  reaction  for  which  there  is  a  net  increase  in  particles  has  an  equilib- 
rium constant  which  decreases  with  density.     An  example  is 

0  =  ~2  °2 

where  two  oxygen  atoms  are  produced  for  each  molecule  of  molecular  oxygen.  The 
net  result  of  this  combination  of  enhanced  effective  equilibrium  constant  with 
density  for  reactions  in  which  the  number  of  particles  is  decreased  and  decreas- 
ed effective  equilibrium  constant  with  density  for  those  in  which  it  is  increas- 
ed is  to  produce  a  decrease  in  the  total  number  of  particles  with  density  -  a 
result  which  is  clearly  visible  in  Table  3. 

With  very  few  exceptions  (e.g.  some  of  the  ionic  species)  all  virial 
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coefficients  are  expected  to  be  positive  at  the  temperatures  of  interest  (and, 
in  fact,  even  at  much  lower  temperatures).     It  is  for  this  reason  that  the  HS 
model  is  expected  to  be  particularly  appropriate  here.     Because  of  this  positive 
behavior,  at  these  temperatures  density  corrections  to  the  equation  of  state 
must  be  positive  at  all  densities  and  are  required  to  increase  with  density. 
Any  approximation  which  has  this  character  must  therefore  produce  a  density 
correction  to  the  ideal  gas  which  has  the  correct  sign,  if  not  the  proper 
magnitude.     The  simplest  such  correction  is  associated  with  taking  only  the 
second  virial  coefficient.     According  to  Table  3,  such  an  approximation  does 
produce  an  increase  in  compressibility  factor  with  density  for  the  intermolecular 
potential  functions  used  here  since  these  produce  positive  second  virial 
coefficients.     The  decrease  in  the  magnitude  of  this  correction  with  increasing 
temperature  at  constant  density  is  caused  by  a  reduction  with  increasing  temp- 
erature in  the  magnitude  of  the  second  virial  coefficient  for  many  of  the  inter- 
actions used  since  the  temperatures  of  interest  are  above  those  at  which  the 
second  virial  coefficients  for  these  interactions  exhibit  maxima. 

The  GB  approximation  is  strongly  dependent  on  our  second  virial  coeffi- 
cients, since  those  authors  made  use  of  our  earlier  results  for  the  second 
virial  gas  to  tie  down  their  extrapolations  at  5000  and  6000  kelvins.     For  this 
reason,  their  predictions  should  be  in  close  agreement  with  the  results  for  the 
second  virial  gas  up  to  the  densities  at  which  that  approximation  is  expected 
to  be  valid  or  to  densities  somewhat  below  100  amagats.     This  is  essentially 
the  behavior  exhibited  in  Table  3.     Since  the  GB  model  is  based  on  the  results 
of  our  second  virial  approximation,  this  agreement  is  not  a  test  of  the  GB  model 
but  rather  serves  as  a  test  of  the  computer  program  developed  by  those  authors 
as  well  as  of  our  own  and  as  a  test  of  the  GB  input  data  as  obtained  from  our 
virials.     Since  their  approximation  includes  an  estimate  for  the  effect  of  third 
virial  coefficients  and  since  the  HS  approximation  does  also,  and  since  both  of 
these  are  of  the  same  sign,  the  predictions  of  the  GB  calculation  should  be 
expected  to  agree  with  results  for  the  HS  model  to  slightly  higher  densities 
than  the  second  virial  gas  (and  this  is  also  exhibited  in  Table  3). 

Our  present  results  are  entirely  compatible  with  the  predictions  of  the 
second  virial  gas.      This  is  no  more  than  expected  since  both  calculations  were 
based  on  the  same  second  virial  coefficients  and  since  our  model  reduces 
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exactly  to  the  second  virial  gas  in  the  limit  of  low  density  when  all  but 
linear  density  terms  are  neglected. 

As  the  density  increases  beyond  approximately  100  amagats,  our  results  for 
the  HS  model  begin  to  deviate  very  rapidly  from  the  second  virial  results, 
becoming  over  twice  as  large  at  a  density  of  1000  times  normal.     It  should  be 
noted  that  there  is  a  similar  relationship  between  the  results  for  the  GB 
approximation  and  those  for  the  second  virial  gas  except  that  since  those  two 
approaches  involve  very  similar  approximations,  differences  between  their 
results  are  considerably  smaller  than  between  the  HS  and  second  virial 
approaches,  especially  at  the  highest  densities. 

C.       SUGGESTED  POSSIBILITIES  FOR  IMPROVEMENTS  IN  THE  TABLES. 

As  mentioned  above,  a  number  of  approximations  were  made  in  order  to 
expedite  the  completion  of  these  calculations.     These  approximations  were 
considered  to  be  sufficiently  minor  so  as  not  to  affect  a  study  of  the  effect 
of  density  corrections  on  the  concentrations  and  of  the  feasibility  of  carrying 
out  such  calculations.     Because  of  the  nature  of  the  calculation  and  because  of 
the  applicability  of  Le  Chatelier's  principle,  most  of  the  approximations  should 
not  be  expected  to  affect  the  accuracy  of  the  compressibility  factors  drastic- 
ally.    In  this  section  we  shall  describe  ways  in  which  these  approximations 
might  be  relaxed  in  order  to  produce  more  accurate  tables.     Please  note  that 
the  order  in  which  the  approximations  are  discussed  is  not  necessarily  related 
to  the  order  of  their  importance. 

In  this  calculation  we  include  only  the  compressibility  factor  from  among 
the  thermodynamic  properties.     Thermodynamic  properties  can  be  obtained  either 
through  numerical  operations  on  our  tables  of  compressibility  factors  or 
through  direct  calculation  of  the  properties  from  the  equations  given  in  the 
text.     For  some  of  the  properties,  such  direct  calculations  require  a  knowledge 
of  the  temperature  derivatives  of  the  temperature  dependent  parameters  associat- 
ed with  the  effective  potential  function,  i.e.  b(T)  and  e(T),  particularly  the 
first.     This,  in  turn,  requires  an  improvement  in  the  method  used  for  the  cal- 
culation of  the  hard  sphere  diameters  for  temperatures  above  that  at  which  the 
second  virial  coefficient  attains  its  maximum.     While  the  method  used  by  us 
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produces  a  relatively  trivial  discontinuity  in  the  value  of  b(T)  at  the  change- 
over temperature,  it  does  produce  a  discontinuity  in  its  temperature  derivative 

which  is  much  larger.     It  is  clearly  possible  that  these  discontinuities  in  the 
db.  . 

values  of      jj     will  also  be  reflected  in  dZ/dT  values  obtained  in  a  numerical 


differentiation.     Hence,  even  though  it  might  be  possible  to  neglect  the  temp- 
erature derivative  of  b(T)  as  being  small,  it  is  probably  very  necessary  to 
improve  on  the  method  of  calculating  b(T)  at  high  temperatures  to  ensure  that 
spurious  discontinuities  in  slope  are  not  introduced. 

There  are  a  number  of  methods  which  could  be  used  for  the  calculation  of 
b(T)   for  temperatures  above  that  at  which  the  second  virial  attains  a  maximum 
and  these  should  be  investigated.     Since  the  effect  of  attraction  can  be 
totally  neglected  at  these  temperatures,  the  most  promising  method  might  be  one 
in  which  the  problem  of  finding  the  two  parameters  b(T)  and  e(T)  at  each  temp- 
erature is  replaced  by  that  of  obtaining  b(T)  only.     The  present  method 
essentially  does  this  but  in  a  very  arbitrary  manner,  and  must  be  modified  so 
as  to  produce  a  smooth  table  of  values  for  db/dT. 

An  obvious  improvement  in  the  tables  will  also  result  when  the  inter- 
molecular  potentials  used  by  us  are  replaced  by  improved  ones.     The  most 
important  of  these  have  already  been  determined  by  us  in  earlier  work  under 
this  contract.     In  an  earlier  report,  we  estimated  the  possible  effect  of  this 
on  calculated  tables.     Although  the  effect  was  shown  to  be  considerably  smaller 
than  were  the  substantial  differences  reported  here  between  our  results  and  those 
of  the  GB  model,  they  were  nevertheless  found  to  be  not  negligible.     Such  a 
study  needs  to  be  made  within  the  context  of  the  HS  theory,  it  being  otherwise 
impossible  to  place  meaningful  estimates  of  precision  on  our  results.     For  the 
second  virial  coefficient  alone,   the  ratio  of  the  value  predicted  for  the  (18,6) 
potential  to  that  predicted  for  the  (12,6)  at  a  temperature  of  5000°  when  the 
value  at  500°  is  correctly  predicted  by  both,   is  given  approximately  by 


dT 


/"5000  V 
V500  / 


2 
18 


12 


=  (10) 


1 
18 
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for  a  difference  of  approximately  14%. 

Those  properties  which  depend  on  db/dT  would  probably  only  be  modified 

slightly  by  changing  the  potentials  since  db/dT  is  small  for  almost  all 

27 

potentials  at  the  temperatures  of  interest.     In  our  earlier  report      this  was 
shown  to  be  true  for  the  first  density  corrections  to  such  properties. 

Improvement  of  the  tables  through  the  improvement  of  certain  of  the  inter- 
molecular  potentials  used  (e.g.  pair  interactions  involving  NC^)  would  take 
considerable  additional  effort.     This  would  require  a  literature  search  for 
experimental  data  for  second  virial  coefficients  and  viscosity  data  for  the 
relevant  species  and  the  determination  of  parameters  for  intermolecular  forces 
using  such  data. 

Related  to  this  but  somewhat  broader  in  scope  is. >the  need  for  a  detailed 
study  of  mixing  rules  by  means  of  which  potential  functions  which  describe  the 
interactions  between  unlike  species  are  inferred  from  those  which  describe  the 
interactions  between  like  species.     This  is  particularly  important  for  such 
pairs  as         NO  which  are  major  constituents  under  the  conditions  of  interest. 

One  of  the  unexpected  problems  which  we  met  in  this  work  had  to  do  with 
the  need  for  having  enough  interaction  virials  for  the  description  of  the  inter- 
actions between  a  given  major  species  and  other  important  species.     Because  the 
net  effect  of  these  interactions  on  the  equilibrium  constant  generally  appears 
as  a  smaller  difference  between  larger  quantities,  such  quantities  cannot  be 
arbitrarily  neglected.    We  solved  this  problem  partly  by  a  shift  to  other 
reference  species  and  partly  by  arbitrary  approximation  of  the  unknown  inter- 
actions.    A  study  needs  to  be  made  to  establish  a  criterion  for  determining 
when  such  interactions  can  be  neglected.     The  importance  of  this  can  be  seen 
from  the  sheer  number  of  possible  interactions  which  can  be  needed  in  a  cal- 
culation of  this  kind.     Thus,  in  a  mixture  of  n  constituents  there  are 
n(n+l)/2  pair  interactions.     A  mixture  of  30  species  therefore  has  465  possible 
pair  interactions!!     The  determination  of  these  would  constitute  a  tremendous 
job  made  particularly  difficult  in  our  case  by  the  facts  that  data  are  not 
available  for  the  appropriate  binary  mixtures  (since  many  of  these  species  can- 
not be  handled  at  ordinary  temperatures  while  other  species  are  not  available 
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in  more  than  trivial  concentrations  at  equilibrium  at  such  temperatures).  This 
is,  for  example,  true  of  the  atomic  species.     Another  major  problem  comes  from 
the  fact  that  of  the  465  pairs,  435  involve  interactions  between  unlike  species. 
The  possibility  of  there  being  data  on  the  435  binary  mixtures  needed  at  a 
sufficient  number  of  temperatures  and  for  a  sufficient  number  of  relative 
concentrations  from  which  to  infer  potentials  of  interaction  is  extremely  small. 
Thus,  many  of  the  interactions  need  to  be  estimated  by  whatever  means  is 
available.     Clearly,  any  reduction  in  the  number  of  pair  interactions  needed 
and  in  the  accuracy  with  which  the  remaining  ones  are  needed  produces  a 
comparable  direct  reduction  in  the  amount  of  work  involved  in  calculations  of 
this  kind.     Thus,  criteria  need  to  be  established,  within  the  framework  of  our 
model,  by  means  of  which  it  can  be  determined  when  a  particular  pair  inter- 
action can  be  neglected  and  when  a  pair  interaction  contributes  sufficiently 
little  so  that  it  can  be  approximated  in  a  rather  cavalier  fashion. 
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Fig,  2,  Effective  potentials  vs  molecular  separation  in  arbitrary  units  for  several 

VALUES  OF  M  AND  FOR  T*  =  2,5,     It  IS  SEEN  THAT  THOUGH  THE  SHAPE  OF  THE  "BOWL" 
IS  QUITE  SENSITIVE  TO  THE  CHARACTERISTIC  PARAMETER  M,  THE  SPHERE  DIAMETER  AT 
THIS  TEMPERATURE  IS  ONLY  WEAKLY  SENSITIVE  TO  IT, 
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.2        .4        .6  .8 

p* 


Fig,  4.  Compressibility  factor  for  argon  vs  reduced  density  pGV)  =  p*  for 
T  =  119. SIC.  Curves  L  this  work;  Curves  2  and  l\,  ZBH  theory  from 
references  22  and  23,  respectively;  Curves  3,  ZBH  theory  based  on  Monte 
(Precalculations  of  perturbation  terms     Curves  5,  FVT  experimental 
datajj.  The  "computer  experiments"^  are  designated  by  circles,  Monte 
Carlo;  and  squares,  molecular  dynamics. 
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Fig,  7,  Compressibility  factor  vs  density  for 

ARGON  FOR  T  =  328 , 25Kj  SEE  CAPTION/  FlG .  4, 


65 


AEDC-TR-76-85 


AEDC-TR-76-85 


5.0 


ARGON 
T=  673K 


1.0 


.4 


.6 


.8 


Fig,  9,  Compressibility  factor  vs  density  for 
argon  for  T  =  673K;  see  caption,  Fig,  4, 
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Fig,  10. 


Sensitivity  of  equation  of  state  for  3sH  theory  to  choice  of  potential 
function.  Compressibility  factor  is  plotted  vs  density  in  amasat.  For 
Curve  1,  the  reference  potential  is  the  Lennard-Jones  (18,6);  for  Curve  2, 
the  Lennard-Jones  (12,0 ,  The  numerical  results  are  from  calculations 

FURNISHED  BY  TOXVAERIT  USING  POTENTIAL  PARAMETERS  e/K  =  119.8,  °  =  3.VB 
FOR  THE  LATTEF  AND  EQUATION  (12)  FOR  THE  FORMER. 
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Fig.  11,    Sensitivity  of  present  theory  to  choice  of  potential  functions.  Potential 
functions  are  as  described  in  caption  for  flg.  10.  the  small  difference 
between  Curves  1  and  2  is  mostly  due  to  the  slight  differences  in  second 

VI RIALS  PRODUCED  BY  THE  TWO  POTENTIALS  IN  THE  VICINITY  OF  THE  ISOTHERM. 
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z 


yo(Amagat) 

Fig,  12,  Compressibility  factor  vs  density  in 
amagat  for  nitrogen  t  *  123, 15k,  curves  1/ 
present  theorxi  curves  2,  experimental  pvt 
measurements. 
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Fig.  13,  Compressibility  factor  vs  density  in 
amagat  for  nitrogen  t  =  173,1510  see  caption/ 
Fig.  12. 
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Fie.  14,  Compressibility  factor  vs  density  in 
amagat  for  nitrogen  t  =  273. 1510  see  caption, 
Fig.  12. 
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Fig,  15,  Compressibility  factor  vs  density  in 
^magat  for  nitrogen  t  =  675k;  see  caption, 

ZIG,  12, 
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Fig,  16,  Concentration  versus  density  for 
several  important  species  for  t  =  3000k, 
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Fig,  17,  Concentration  versus  density  for 
several  important  species  for  t  =  900ok, 
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